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ABSTRACT 


An approach for computing worst-case flutter margins has been formulated in a robust stability 
framework. Uncertainty operators are included with a linear model to describe modeling errors and flight 
variations. The structured singular value, p, computes a stability margin that directly accounts for these 
uncertainties. This approach introduces a new method of computing flutter margins and an associated 
new parameter for describing these margins. The (I margins are robust margins that indicate worst-case 
stability estimates with respect to the defined uncertainty. Worst-case flutter margins are computed for 
the F/A-18 Systems Research Aircraft using uncertainty sets generated by flight data analysis. The robust 
margins demonstrate flight conditions for flutter may lie closer to the flight envelope than previously 
estimated by p-k analysis. 
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CHAPTER 1 
INTRODUCTION 


Aeroelastic flutter is a potentially destructive instability resulting from an interaction between aero- 
dynamic, inertial, and structural forces. 1 The stability properties of the aeroelastic dynamics must be 
investigated to determine a flight envelope that is clear of flutter instabilities for new aircraft designs or 
new configurations of current aircraft. Analytical predictions of the onset of flutter must be accurate to 
reduce dangers and costs associated with experimental estimation. 

Critical flutter conditions are the points closest to the flight envelope at which flutter instabilities 
occur. This concept of closeness is formally defined here as a flutter pressure that considers the critical 
dynamic pressure for a constant Mach value. Obviously, different flutter measures such as a flutter speed 
can be defined because a unique equivalent airspeed is associated with each dynamic pressure for a given 
Mach number; however, this paper uses definition 1.0.1 for a flutter pressure to describe the critical 

flutter flight conditions. 

Definition 1.0.1: Aflutter pressure is the smallest value of dynamic pressure for which an aircraft at 
a particular Mach number experiences aflutter instability. 

The flutter pressure is used to compute a stability margin, or flutter margin, that indicates the distance 
between the flutter pressure and a reference point. A common flutter margin, T, considers the difference 
in dynamic pressure between the flutter pressure and a point on the edge of the flight envelope. Another 
common flutter margin, n, considers the percentage difference between equivalent airspeeds at the flutter 
condition and a point within the flight envelope. 

Definition 1.0.2: Aflutter margin relates a measure of distance between the flight condition associ- 
ated with the flutter pressure and a reference point. 

The traditional p-k method has been extensively used to compute flutter margins for a variety of 
military and commercial aircraft. 5 This iterative method uses an analytical dynamic model coupled with 
harmonic motion solutions for the unsteady aerodynamic forces. The p-k method predicts flutter margins 
entirely from a theoretical model that may not accurately describe the true dynamics of the airplane. The 
resulting flutter margins do not account for possible variations between the model and the aircraft. 

The community studying aeroelasticity has identified the development of improved methods for 
characterizing flutter margins as a vital research area. 4 Flight testing for envelope expansion incurs 
dramatic time and cost because stability margins are not computed with a high level of confidence using 
traditional methods. 2 The flutter dynamics often exhibit an explosive behavior that results in a sudden 
change in stability for a small change in flight conditions. Thus, small errors in predicted margins could 
have grave consequences for aircraft and crews operating near the flutter conditions. 

Several approaches exist for characterizing accurate flutter margins using flight data generated by the 
aircraft. These data describe the true dynamics and can be used to generate realistic models and compute 
confident flutter margins. Parameter estimation algorithms have been developed to directly identify an 
aeroelastic model from the flight data. 5 ’ 6 The accuracy of the resulting model can deteriorate as the com- 
plexity and number of degrees of freedom of the system increase and signal-to-noise ratios decrease from 
optimal wind-tunnel conditions to realistic flight levels. Modal filtering has been introduced in associa- 
tion with parameter estimation algorithms to simplify analysis by decoupling the system into a set of 
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first-order responses. ’ This type of filtering does not guarantee robustness and may not perform well 
for systems with many closely-spaced modal natural frequencies that cross and shift as flight 
conditions change. 

Other approaches toward computing confident flutter margins evaluate the robustness of a stability 
margin with respect to changes in the model as an indication of the confidence in that margin. A flutter 
margin robust to perturbations to the model is a confident margin because model inaccuracies do not 
affect that margin. An algorithm has been developed to compute the most critical flutter margin with 
respect to first-order perturbations in a model. 9 This method considers only parametric perturbations and 
can be computationally expensive. A robust control framework has been adopted using a feedback struc- 
ture to relate the structural model and the aerodynamic model. 10 This approach uses highly conservative 
robustness conditions with respect to an uncertainty structure that may not be physically meaningful. A 
similar approach is adopted allowing unmodeled dynamics and high-order parametric perturbations 
based on series expansion. 11 Statistical approaches are also considered to formulate a flutter probability 
measure. 12, 13 These approaches will converge to a robustness indicator using Monte Carlo simulations, 
but the computation time can be prohibitive for complex systems. The robustness measures for these 
perturbation and statistical approaches are suspect because no global guarantees can be made as to pertur- 
bations not explicitly considered by the minimization algorithms or the Monte Carlo simulations. 

An approach to computing flutter margins that guarantees a level of robustness and directly accounts 
for flight data is presented herein. 14 An aeroelastic model is formulated in a formal robust stability frame- 
work that uses a set of norm-bounded operators. A, to describe modeling errors and uncertainty. A multi- 
variable robust stability measure known as the structured singular value, |i, computes flutter pressures 
that are robust to the amount of modeling errors as determined by A . 15 A robust flutter margin problem is 
posed by questioning what is the largest increase in dynamic pressure for which the plant is stable despite 
possible modeling errors described by A. 

Flight data are easily incorporated into the ji analysis procedure. The modeling errors are determined 
by comparing transfer functions obtained by flight data with transfer functions predicted by the analytical 
model. The norm bound on A is chosen based on these observed errors. A model validation condition is 
used to ensure the A is sufficient to account for multiple data sets without being excessively conservative. 
With respect to A, a worst-case flutter boundary is computed that directly accounts for flight data. 

This method is inherently different from traditional algorithms based on p-k methods or parameter 
identification and robustness approaches. The (I method uses information from both an analytical model 
and flight data; traditional approaches use only one of these sources (fig. 1.1). Methods that use only an 
analytical model can be inaccurate, and methods that use only the flight data can fail if the data are of 
poor quality. The (I method uses the flight data to improve the analytical model by adding uncertainty 
operators. Poor quality flight data will merely increase the difficulty of obtaining a reasonable uncertainty 
description resulting in a small A. The robust margins will be similar to the nominal margins in this case, 
which makes intuitive sense because any information obtained from the data should only enhance the 
plant model and improve the accuracy of the flutter margin. 

The concept of computing robustness in flutter margins has been recognized for its importance and 
has recently been termed a state-of-the art research area in aeroelasticity. 16 Informal measures of robust- 
ness are not necessarily useful because the informal measures provide no guarantee for the system stabil- 
ity. The p method is based on operator theory and provides a well-defined concept of robustness that has 
a clear set of guarantees as to the stability properties of the system. 
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Figure 1.1. Information flowchart for traditional and p methods to compute flutter margins. 


The p method also allows consideration of a rich set of modeling errors with direct physical interpre- 
tations. Traditional techniques compute robustness with respect to a theoretical estimate of modeling 
errors that can be too conservative or overly optimistic. The jx method is able to use flight data to deter- 
mine these errors by employing model validation algorithms that are developed specifically for compar- 
ing data with uncertain models in the p framework. Thus, the resulting stability margins directly account 
for the flight data when determining robust flutter pressures. 

Robust, or worst-case, flutter margins are computed for the F/A-18 Systems Research Aircraft (SRA) 
using the p method. The flutter results presented represent an improvement to previous flutter results for 
the F/A-18 SRA computed using the traditional p-k method. Nominal flutter margins computed using the 
p method but ignoring all uncertainty operators are shown to closely match the p-k method flutter 
margins. These results lend validity to the p method as an accurate indicator of flutter instability. Robust 
flutter margins are computed with the p method by directly accounting for modeling uncertainty to repre- 
sent conservative worst-case stability margins. The computation of the p method takes only a few min- 
utes on a standard workstation. 
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CHAPTER 2 
ROBUST STABILITY 


2.1 Signals 

A common definition of a signal is a Lebesgue measurable function that maps the space of real 
numbers R into R n . A space of such signals is denoted S. 

Definition 2.1.1: The space of signals that are Lebesgue measurable functions is S- 

S = {f: R->R n } (1) 

Analog measurements x(t) of physical systems are real vector functions of the real parameter 
t describing time and thus are valid members of the space of signals, x(t) e S. Values of the time param- 
eter, which are often arbitrarily numbered as a distance from some reference point, actually extend to 
positive and negative infinity. Stability for physical systems must ensure stability for all values of time. A 
time-domain 2-norm is defined as a measure of size (or energy) for time-domain signals that considers 
all time. 

Definition 2.1.2: The 2-norm measures the energy of a signal x(t) e S. 

I 

H x(t )|| 2 = (CJ x(t) ^ dt T (2) 

One characteristic of a stable system is only finite-energy output signals are generated in response to 
finite-energy input signals. Signals with finite energy are known as “square integrable” because the inte- 
gral of the square of the signal is finite. The Lebesgue space of square integrable signals is defined as 
^ 2 (-°°,°°). This space is also referred to as the infinite-horizon Lebesgue 2-space to denote that the norm 
uses an integral over infinite time. 

Definition 2.1.3: The space ^ 2 (— ° 0 , 00 ) consists of square integrable time-domain signals. 

^ 2 (-°°,°°) = (x(t) :xe5, ||x(t)|| 2 < oo} ( 3 ) 

Signals associated with physical systems are only known for values of time greater than the time at 
which measurements are started. Stability analysis and norm computations using these signals cannot use 
properties of the signal before the starting time because no information is known. The traditional method 
of characterizing these signals is to assume the signal is identically zero for all times before the starting 
time. The time value at which measurements are started can be chosen without loss of generality and is 
usually chosen to be t = 0. The space ^ 2 [0>°°) is defined as a subset of A(-oo oo) to emphasize 
such signals. 2 

Definition 2.1.4: The space ^ 2 [0,°o) c ^ 2 (-°°,°°) consists of square integrable time-domain 
signals that are identically zero for all t < 0. 

^ 2 [°^) = |x(t) : j^[|x(t)|| 2 dt<oo, x(t) = 0 for all t < 0 
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A similar space Z,(-~ 01 c uf,(- is defined for signals x(t) that are assumed to begin at 
t = _oo and are identically zero for all times t = 0. The integral to compute the energy for elements of 

this space considers t = until t > 0. 

Frequency-domain signals are often considered in stability analysis but do not fallinto the set of 
signals, S. These signals/C; CO) are complex-valued functions of the imaginary untty = J- 1 , and the real 
frequency variable o> is expressed in rad/sec. The set, S } „ , is defined for frequency-domatn s, goals. 

Definition 2.1.5: The space of frequency-domain signals is S jQ) . 


Sj a = {/(;' CO) : j R C n and /*0'co) = /H®)} ^ 

A frequency-domain 2-norm is formulated to compute a measure of energy. 

Definition 2.1.6: The 2-norm measures the energy of the signal f(j (0) e Sj or 

1 

l/lj = (jj fj’UWUWef («) 

A frequency-domain Lebesgue space, -d 2 , is defined for finite-energy signals. 

Definition 2.1.7: The space J 2 consists of frequency-domain signals with finite energy. 

= {/O’©) : / G Sjvi' H/H2 < 00 ^ (7) 

The spaces ^ 9 (-oo,oo) and are isomorphic Hilbert spaces under the appropriate inner products 
through the Fourier transform, which means the spaces have equivalent algebraic properties. This 
relationship is used to simplify notation by rarely distinguishing between time-domain and frequency- 
domain signals except where the context does not make it clear. The notations for the 2-norm of time 
domain and frequency domain are also not distinguished because the notations are equivalent, as demon- 
strated by Parseval’s identity. 17 

An important subspace of J 2 is the Hardy space, » 2 . This space contains the complex variable func- 
tions that are analytic in the open right-half of the complex plane and have finite 2-norm. 

Definition 2.1.8: The Hardy space, c J 2 , consists of the following functions. 

‘» 2 = {/( s) : /( s) G J 2 and f (s) is analytic in Re(s) >0} ( 8 ) 

The space A[0,oo) is isomorphic to ^ 2 through the Laplace transform. Members of this set can be 
interpreted as frequency-domain representations of finite-energy time-domain signals that are assumed to 

be zero at t < 0. 


2.2 Systems 

A system P is defined as an operator mapping the space of input signals S in to the space of output 
signals S our This definition implies that for any w e S in and z = Pw, then z € S our 
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( 9 ) 


P • C. — > £ 
r ' ^in J out 

This paper will consider linear, time-invariant systems defined by state-space equations. 

x(t) 

_y(t)_ 

The signal x(t) e R ns is the state vector, u(t) e R n ‘ is the input vector, and y(t) e R n ° is the output 
vector. The state update matrix is A p e R^*" 5 ; B p e R n s Xn i determines how the input affects the 
states; C p e R n o Xn s computes the outputs as a linear combination of states; and D p e R 0 1 is the 

direct feedthrough from inputs to outputs. The operator S = {A p , B p , C p , D p } denotes the time-domain 
state system. 

Linear time-invariant state-space systems are commonly represented by transfer-function operators. 
These functions, P(s), are complex-valued matrices of the complex Laplace transform variable, s. Such a 
transfer-function matrix exists if and only if the state-space system is linear and time-invariant. 

P(s) = Dp + Cp(sI ns -A p )- 1 Bp (11) 

Stability must be considered over the infinite-horizon time lengths so that the operators used map 
.^ 2 (-oo,oo) into ^ 2 (-co,oo). Properties of the Fourier transform relating ^t -00 ’ 00 ) an d ^2 i m piy a state- 
space system, S : ^f -00 ’ 00 ) ^li- 00 , 00 ), is linear and time-invariant if and only if the associated 

transfer- function matrix P is such that y = Pu e f° r any u e This condition leads to consideration 
of the gain for these signals. 


A p B p 
C p D p 


x(t) 

u(t) 


( 10 ) 


M2 m 2 

I|u|| 2 " I|u|| 2 u ; 

This ratio of 2-norms will be finite if the system is stable. Properties of the 2-norm are used to derive 
a condition on the system transfer-function operator, P. This condition is referred to as an induced norm 
because it results from consideration of signal norms associated with the operator. The ^f 2 induced norm 
is defined as the # -norm. 

oo 

Definition 2.2.1: Define the ^-norm for transfer-function operators. 

I|P|L = sup G(P(;' 0 )) (13) 

0 

A space of operators with finite 3^-norm is denoted as . 

Definition 2.2.2: The space consists of systems with finite ‘& oo -norm. 


= {P: «PIL<~} (14) 

Transfer functions of linear time-invariant systems are stable if and only if z = Pw and w e implies 
z € ^ 2 . This implication results from the Laplace isomorphism between ^ 2 [0,°o) and # 2 space. These 
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transfer functions are shown to be analytic in the open right-half complex plane with finite norm. 
Define the space to contain these operators. 

Definition 2.2.3: The space consists of transfer Junctions of stable, linear, time-invariant 
systems with finite ‘& oo -norm. 

= {P : P is analytic in Re(s) > 0 and llPlL < °°} ( 15 ) 

A subspace is often defined for rational elements. 

Definition 2.2.4: The space, c ^ , consists of rational elements of ^ . 

= {P : P e ^ and P is rational } (16) 

Transfer-function operators of linear, time-invariant state-space systems are rational functions of the 
Laplace transform variable, s. These transfer functions P e if and only if P is stable such that no 

poles lie in the closed right-half plane. The space which may appear to be a mathematical abstrac- 

tion, is thus shown to have a physical interpretation. is merely the operator theory representation of 
stable, rational, transfer functions. 

2.3 Small Gain Theorem 

Stability of a linear time-invariant system is determined by location of all poles in the left-half plane. 
Robust stability in the and p frameworks is determined by considering an interconnection of stable 
operators. The basis for” determining stability of these interconnections of operators is the “small 
gain theorem.” 

The small gain theorem states that a closed-loop feedback system of stable operators is internally 
stable if the loop gain of those operators is stable and bounded by unity. Several formulations of the small 
gain theorem are derived for various signals and systems. Theorem 2.3.1 presents the formulation used 

for this paper. 17 

Theorem 2.3.1 (Small Gain Theorem): Given the feedback interconnection structure of figure 2.1 
for stable transfer-function operators P, A : ^ G if W^-norm of the loop gain 

is bounded by unity such that flPAH^ < 1 , then: 

1. the closed-loop system is well-posed and internally stable. 

2. a unique y, w e A 2 I s associated with each u e • 

This small gain theorem is overly restrictive in the sense of requiring P, A e A more general 

small gain theorem is formulated for operators not restricted to lie in the subspace theorem 2.3.1 is 

a special case of this general theorem. 18 The extended operator space in the general small gain theorem 
allows consideration of robustness for systems composed of nonlinear and time-varying operators. The 
requirement of considering stable, rational, transfer-function operators is explicitly stated in the theorem 
to emphasize that the nominal aeroelastic system considered in this paper is assumed to be stable and the 
flutter margin is associated with a destabilizing perturbation to that nominal system. 
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Figure 2. 1 . Block diagram for the small gain theorem. 


The second condition in theorem 2.3.1 is associated with the first condition guaranteeing a well-posed 
and stable system. This uniqueness condition can be understood by consideration of the solution y for the 
loop equations shown in figure 2.1. 


y = P(u + w) 

-1 07) 

= (I - PA) Pu 

The inverse term, (I - PA) , has a magnitude of infinity if the norm of PA is allowed to be unity. 
Such a condition would allow the norm of signal y to be infinite despite a norm-bounded u input signal. 
Restricting ||P Al)^ < 1 ensures the inverse term exists and a unique finite-norm y is generated in response 
to a finite-norm u. The issue of well-posedness requires this condition to hold at s = °o and is automati- 
cally considered by the -norm. 

2.4 Robust Stability 


The small gain theorem can be directly used to analyze robust stability of a plant model with respect 
to a set of perturbations. These perturbations are used to describe uncertainty in the analytical plant 
model caused by errors and unmodeled dynamics. Usually, the exact value of the modeling error is not 
known, but a norm-bounded, real scalar, (X > 0, can be placed on the size of that error. Define the A 
of norm-bounded operators describing these perturbations that affect the plant P through a 
feedback relationship. ° 


A= (A : || A|| ^ < a} (18) 

The small gain theorem allows consideration of the entire set of possible modeling uncertainties as 
described by all A e A. The -norm of the loop gain cannot be explicitly computed for these systems 
because an infinite number of loop gains P A generated by the A exists. The triangle inequality of norms 
can be used to generate a sufficient condition for robust stability of P. 

IIPAIL < IIPILIIAIL (19) 


A condition for robust stability of the closed-loop system can be stated. 

Lemma 2.4.1: The plant P is robustly stable to the set of uncertainty perturbations, A that enter the 
system as in figure 2.1 with [| A|| ^ < a for all A e A if 
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( 20 ) 


«PlL < S 

I emma 2 4 1 shows a sufficient, but not necessary, condition for robust stability. Th e structured 
singular value. |X, is iniroduced in the nett chapter as a less conservative measure of robust stabthty that 

sufficient and necessary. 

An excellent illustrative example has previously been presented 18 to demonstrate the issue of robust 

s^enUhat^is S seento > be WenticaHo tte g^rated^sing the small gain theorem 

Sd Tmma 2.4.1. A similar example is given below for the feedback interconnection in figure 2.2. 
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Figure 2.2. Block diagram with uncertainty for the example system. 


The single-input and single-output elements in the nominal system model (fig. 2.2) are p, whic 
represents the plant dynamics; a, which represents actuator dynamics; and k, which represents a feed bac 
controller Each of the nominal system elements are stable transfer functions contained in . A 
modding emor exists on the output of the actuator that is represented by a multiplicative uncertainty 

operator, 5 e , on the output of the element a. 


The transfer function from w to z can be computed as follows. 

akp) 'akp^jw 


z = -(1 + 


( 21 ) 


Internal stability of the closed-loop feedback system is equivalent to stability of the feedback system 
shown in figure 2.3 with the operator g = -( 1 + akp) akp . 



Z 


w 
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Figure 2.3. Block diagram for robust stability analysis of the example system using the small gain 


theorem. 
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Because the operators 8, g <= are stable, the Nyquist criterion determines the closed-loop 

system is stable if and only if the Nyquist plot of 5g does not encircle the -1 point. This stability condi- 
tion is equivalent to the following norm condition. 

sup IgO'co) 5(y'co)| < 1 (22) 

This condition is an fl^-norm condition on the loop gain, g8. Thus, classical Nyquist arguments 
derive an r^-norm condition that is equivalent to the stability condition immediately formulated by 
applying the small gain theorem. 

closed-loop stability Ilg8|| 00 <l (23) 

The error in the actuator command is unknown and possibly time-varying, so the operator 5 is used to 
allow consideration of a range of errors. Assume the actuator is weighted such that the range of errors is 
described by the set of perturbations, 11811^ < 1. Lemma 2.4.1 is used to generate a condition that ensures 
the system is robustly stable to all actuators errors described by 8. 

closed-loop stability <= ||g|| 00 <l (24) 
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CHAPTER 3 

STRUCTURED SINGULAR VALUE 


3.1 Linear Fractional Transformations 

The linear fractional transformation (LFT) is a common framework suitable for robust stability 
analysis using arguments based on the small gain theorem. An LFT is an interconnection of operators 
arranged in a feedback configuration. These operators can be constant matrices, time-domam state -space 
systems, or frequency-varying transfer functions. Consider a linear operator P e C that 

is partitioned into four elements. 


P = 


P 11 P 12 
P 21 P 22 


(25) 


The LFT, F (P, A) , is defined as the interconnection matrix such that the upper loop of P is closed 


with the operator A e C 


'i xo i 


Definition 3.1.1: Given P e C 
loop LFT of? closed with A such that y = F U (P, A) u as in figure 3.1. 


(o,+o 2 )x(i l + i 2 ) an(J Ae q {X 0 \ define F U (P, A )) as the upper- 


- 1 , 


F U (P,A) = P 22 + P 21 A(I-P n A) P 12 


(26) 
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Figure 3.1 Linear fractional transformation F U (P, A). 


A similar LFT is defined as Fj(P, A) to represent the interconnection matrix of the lower loop of 

i 2 x o 2 

P closed with an operator A e C 


(o, +o 2 )x(i, +i 2 ) r ,h xo 2 

Definition 3.1.2: Given Pe C and Ae C 

LFT of? closed with A such that y = F,(P, A) u as in figure 3.2. 


, define F U (P, A) as the lower-loop 


F,(P, A) = P u +P 12 A(I-P 22 A)' 1 P 21 


(27) 
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Figure 3.2. Linear fractional transformation F|(P, A). 


An example of an interconnection that is common in stability analysis is the representation of a time- 
dependent state-space system as a frequency- varying transfer function. Define S as the constant matrix 
whose entries are the {Ap, B p , C p , D p } matrices of a state-space realization. The transfer function can be 
written as an upper-loop LFT involving S and the Laplace transform variable s where | is an integrator. 


x = A p x + B u A p B p 

p p => S = p p 
y = C p x + D p u C p D p 


P(s) = D p + C p (sl - A p ) ! B p 

= f u ( s 4) 


(28) 


The LFT is a useful framework for analyzing complex systems with many feedback and series 
interconnections of operators. Property 3.1.3 shows the main property of LFTs that will be used herein. 
This property allows complex systems of several interconnected LFTs to be expressed as an equivalent 
single LFT. The operators of the new LFT are block-structured with blocks composed of the individual 
operators of the LFTs from the original system. 

Property 3.1.3: Feedback and series interconnections of LFTs can be formulated as a single LFT. 

This issue of stability for LFT systems is associated with the concept of a well-posed interconnection. 
Stability analysis based on the small gain theorem given in theorem 2.3.1 can be used to guarantee the 
LFT is stable and well-posed. 


3.2 Structured Uncertainty 


The concept of uncertainty is formulated as a set of norm-bounded operators, A, associated with a 
nominal plant, P, through an LFT feedback relationship. A family of plants, P, arises through consider- 
ation of F U (P, A) for every A e A. The true plant model is assumed to lie within this family of plants. 

Modeling the uncertainty as a norm-bounded operator can lead to overly conservative models. The 
uncertainty description can be made more accurate by including frequency information. Formulating a 
model of a physical system that is accurate at low frequencies but less accurate for representing the 
system response at high frequencies is often possible. A frequency-varying transfer function, W, is gener- 
ally associated with each uncertainty element to describe magnitude and phase uncertainty as it varies 
with frequency. 

Uncertainty can enter a system model in a linear fractional manner in several general ways. Two 
typical types of uncertainty are termed “multiplicative” and “additive” uncertainty. Multiplicative 
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uncertainty can be either on the input or output of a system. Systems with these types of uncertainty are 
easily described in block diagram form. Figure 3.3 shows the LFT for a plant with input multiplicative 
uncertainty. Figure 3.4 shows the plant with output multiplicative uncertainty. Figure 3.5 shows 

additive uncertainty. 



Figure 3.3. Family of plants P = P(I + A W) with input multiplicative uncertainty. 
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Figure 3.4. Family of Plants P = (I + A W)P with output multiplicative uncertainty. 



Figure 3.5. Family of plants P = P + AW with additive uncertainty. 


Uncertainty can also be associated with specific elements of the system. These parametric uncertain- 
ties are usually associated with a system operator in a feedback relationship. The number of input and 
output signals of the system operator is increased to account for the additional feedback signals associ- 
ated with the uncertainty operator. This operation can be demonstrated by considering P generated by a 
system with an unknown pole. 




1 

(s + l)(s + x) 


: x e [2,3] 


(29) 


A norm-bounded, real, scalar, uncertainty parameter 8 can be introduced to account for the possible 
variation in pole value. The set of plants can be written in the LFT framework using this uncertainty oper- 
ator and definition 3.1.1. 
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-1 

-i 


p= • 

F U (P, 5) : P = 

s + 2.5 

1 

s + 2.5 
1 

,||5|L<l,5eR> 



_(s + l)(s + 2.5) (s + l)(s + 2.5)_ 

j 


A complex system with several uncertainty operators can be expressed as an LFT with a single uncer- 
tainty operator using property 3.1.3. This operator is structured as a block-diagonal operator with each 
block associated with the individual uncertainty operators. Two main types of uncertainty blocks exist. A 
full-block uncertainty is a matrix with unknown elements. This type of block is used to describe unstruc- 
tured uncertainty in a group of signals. 

Definition 3.2.1: A full-block uncertainty, A e C n X m , has unknown elements A- for every i e [ 1 , n] 
and] e [1, m], 

A repeated-scalar block introduces more structure into the uncertainty description than a full block 
does. Only the diagonal elements of the matrix contain unknown elements; the remaining elements are 
zero. Furthermore, the diagonal elements are identical. This type of uncertainty is used to relate input- 
output signal pairs with the same uncertainty parameter. 

Definition 3.2.2: A repeated-scalar block uncertainty A e C" X n has zero-valued elements except 
an unknown parameter 8 along the diagonal such that A = 5l n . A scalar block is a repeated- scalar 
block of dimension 1. 

The single structured uncertainty block used for robust stability analysis is formally defined in terms 
of these blocks. Let integers m, n, and p define the number of real scalar, complex scalar, and complex 
full blocks respectively. Define integers Rj, ..., R m such that the i'^ repeated-scalar block of real, para- 
metric uncertainty is of dimension R ; by R ; . Define similar integers Cj, ..., C n to denote the dimension 
of the complex repeated-scalar blocks. The structured uncertainty description A is assumed to be norm- 
bounded and belonging to the following set. 


A = 


(A dia § (X ^1 JCj’ ^C n * A l’ 


.R 


C: XC: 


: 8j e R, 8j e C, Aj € C 1 ' 


(31) 


Real parametric uncertainty is allowed to enter the problem as scalar or repeated-scalar blocks. 
Complex uncertainty enters the problem as scalar, repeated-scalar, or full blocks. Complex uncertainty 
parameters allow uncertainty in magnitude and phase to be modeled; uncertainty in physical characteris- 
tics can be more accurately modeled with real parameters. The robustness analysis will be less 
conservative by accounting for this structure to accurately describe the model uncertainty. 


3.3 Structured Singular Value 

Figure 3.6 shows the general framework for robust stability analysis. The plant operator P(s) e 
is a stable, rational, transfer-function matrix representing the aeroelastic dynamics. A norm-bounded 
A e %& 00 is defined such that A (s) e A describes the modeling errors in P through a feedback relationship. 
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Figure 3.6. Linear fractional transformation system for robust stability analysis using p. 


The robustness of P with respect to the A can be determined using the small gain lhe°rem as 
presented in lemma 2.4.1. This condition guarantees stability for any value A e A if ||P|| • 

robustness condition can be overly conservative because it does not account for structure in the uncer- 
tainty operator. The structured singular value, p, is defined as an alternative measure of robustness. 

Definition 3.3.1: Given the complex transfer-function matrix P e and associated norm- 

bounded set of uncertainty operators A, define p. 


1 

^ (P) " min {a(A) : det(I-PA) =0} 
AeA 


(32) 


Define p = 0 if no A e A exists such that det(I - P A) - 0. 

The structured singular value is an exact measure of robust stability for systems with structured 
uncertainty. The value of p determines the allowable size of uncertainty matrices for which the plant is 


robustly stable as demonstrated in theorem 3.3.2 

Theorem 3.3.2: Given the system in figure 3.6, P is robustly stable with respect to the A, which is 
norm-bounded by real scalar a such that flAll^ < a for alike A if and only if p(P) < - • 

The model P is usually internally weighted such that the range of modeling errors is described by the 
uncertainty set A, which is norm-bounded by 1 . 


|A|| < l for all AeA 


(33) 


Theorem 3.3.3 presents the specific condition for robust stability that will be used in this paper for 
unity norm-bounded uncertainty sets. 

Theorem 3.3.3: Given the system in figure 3.6, P is robustly stable with respect to the A with 
|| A|| ^ < 1 for all AeA if and only i/p(P) < L 

A value of p < 1 implies no perturbation within A exists that will destabilize the feedback system. 
This condition can also be interpreted as saying the true plant dynamics are stable, assuming these 
dynamics lie within the range generated by the nominal model dynamics coupled with the set ot 

modeling errors. 
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Obviously, |i is dependent on the block structure of A. The robust stability properties computed by |i 
will only be accurate if a realistic uncertainty operator is chosen. The structured singular value may be 
arbitrarily greater when computed with respect to an unstructured uncertainty operator as compared to a 
highly structured uncertainty operator. Definition 3.3.1 demonstrates the p condition of theorem 3.3.3 is 
equivalent to the small gain condition of lemma 2.4.1 when the uncertainty is unstructured. 

Unfortunately, p is a difficult quantity to compute. Closed-form solutions exist to exactly compute |l 
for only a small number of block structures for A Upper and lower bounds are used to compute p for 
generalized uncertainty block structures. Appendix A shows a derivation of the upper bound, which rep- 
resents a limit on the worst-case robustness properties. 
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CHAPTER 4 

ROBUST FLUTTER MARGINS 


4.1 Nominal Aeroelastic Model 

21 

Consider the generalized equation of motion for the structural response of the aircraft. 

Mtj + Cti + Kr) + ^Q(s)r) = 0 (34) 

For a system with n modes, define M 6 R as the mass matrix, C 6 R as the damping matrix, 
and K e R n X n as the stiffness matrix. Define q e R as a scalar representing the dynamic pressure, and 
Q(s) e C nX " as the matrix of unsteady aerodynamic forces. This equation is valid for a particular Mach 
number, with a different Q(s) describing the unsteady aerodynamic forces at a different Mach number. 

Values of the aerodynamic force matrix at distinct frequencies can be derived using finite-element 
models of the aircraft and panel methods for unsteady force calculations. This research uses a computer 
program developed for NASA known as STARS. ’ This code solves the subsonic aerodynamic equa- 
tions using the doublet-lattice method. 24 ’ 25 The supersonic forces are generated using a different 
approach known as the constant panel method. 26 

Formulating a linear time-invariant representation of the aerodynamic forces to incorporate them into 
the robust stability framework is desired. Pade approximations can be used to compute a rational function 
approximation to the transfer-function matrix. 


2 S 

Q(s) = Aq + sAj + s A 2 + 


Ao + 


s + (3 1 3 s + P 


(35) 


This form is often referred to as Roger’s form. 27 The equation presented here includes only two lag 
terms, although more terms can be included. The poles of the lag terms, P j and P 2 , are restricted to be 
real and positive to maintain system stability. The matrix elements of Roger’s form can be computed 
using a least-squares algorithm to fit the frequency-varying aerodynamic data. 

The aerodynamic lag terms can be replaced in the formulation with a finite-dimensional state-space 
system represented by a transfer-function matrix using Karpel’s method. 


2 -l 

Q(s) = Aq + sAj+s A2 + Cq(sI — Aq) Bj 


(36) 


Standard system identification algorithms, including curve-fitting or least-squares approaches, can be 
used to compute the elements in the state-space portion of the formulation. The A ( matrices are assumed 
to be known from the low-frequency aerodynamic force data or from experimental wind-tunnel data. 

A matrix fraction approach is also formulated to represent the aerodynamic forces as a linear time- 
invariant system. 29, ^ This generalized form computes rational matrix polynomials in a fractional form 
using a least-squares algorithm. Roger’s form and Karpel’s form can be shown to be subsets of the matrix 
fractional form. 
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The approach taken in this paper is to fit the aerodynamic force matrices to a single, finite- 
dimensional, state-space system. This form is most similar to Karpel’s form, except the additional A- 
matrices are not explicitly accounted for in the formulation. 


Q(s) 


a q Bq 

C Q D QJ 


-l. 


— Dp, + C n (sl Aq) B 


(37) 


Given the number of generalized states, n, and the number of aerodynamic states, n n , define 
n Q x rig n Q x n n x n Q n x n ^ 

A q e R v , Bq e R v , C Q e R v , and D Q e R as the elements of the state-space 
system approximating Q(s). 


Fitting the aerodynamic data to a finite-dimensional state-space system is equivalent to fitting each 
term in the matrix to a real, rational, proper, transfer function. This method seems to contradict the 
methods of Roger and Karpel, which form nonproper transfer functions caused by the terms in s and s . 
An approximation to these forms allows them to fall within the framework of the method used in this 
paper. Including a high-frequency pole in the nonproper term, such as replacing s with - — ^qqqq , would 
not affect the low-frequency region of interest while ensuring stable and proper functions. With the 
approximation, the forms of Roger and Karpel can be shown to be subsets of this method. 

Standard frequency-domain system identification algorithms can generate a system with an arbitrarily 
large number of states. This state dimension does not greatly affect the computational cost of the robust 
stability analysis. Extending the robustness analysis to controller synthesis, however, places an emphasis 
on limiting the state dimension. 51 Limiting the number of states in the identification process is not 
directly considered here, although standard model reduction techniques can be used on the state-space 
system to lower the state dimension. 15 

Generating a state-space representation of the aeroelastic system, including the state-space form of 
the unsteady aerodynarqjc forces, is straightforward. Consider the force vector, y, generated by the state 
vector, T). Define x e R Q as the vector of aerodynamic states. 


y = Q(s)r| o 


X 


a q 

b q 


X 

_y_ 


_ C Q 

d qI 


- 1 !. 


Using x, formulate the aeroelastic differential equation. 


0 = Mt) + Cf) + Kr| + qQ(s)r) 

= Mrj + Cf) + Kr| + qy 
= Mr} + Cf| + Kr| + <?(CqX + DqT() 
= Mfj + Cf) + (K + ^Dq)t| + qCqX 


(38) 


(39) 
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A state-space system is formulated using the generalized states, T| and i], and the unsteady aerody- 
namic states, x. The state-update matrix is determined by the following three differential equations. 


fl 

x 


0 

-M -1 (K + ?Dq) 


-M 'c 

0 


0 


-q M C Q 


11 

X 


(40) 


4.2 Nominal Aeroelastic Model in the Structured Singular Value Framework 

The generalized equation of motion for the nominal aeroelastic system can be expressed in a form 
suitable for using g analysis to compute a flutter margin. The flutter margin is dependent on the flight 
condition parameters that result in a flutter instability, and g is defined to be the smallest perturbation 
among the A that causes an instability. The obvious approach to formulating flutter analysis in the g 
framework is to introduce a perturbation to a flight condition parameter and find the smallest perturbation 

that causes an instability. 

Essentially the two subsystems in the nominal aeroelastic model are composed of the structural 
dynamics, involving mass, damping, and stiffness matrices; and the unsteady aerodynamics scaled by the 
dynamic pressure. The generalized equation of motion demonstrates the dynamic pressure linearly affects 
the dynamics at a constant Mach condition. Perturbations to dynamic pressure can thus enter the system 
through a feedback operator in a linear fractional manner that is perfectly suited to g analysis. 

Consider an additive perturbation, 8- e R, on the nominal dynamic pressure, q Q . 


Separate terms in the system dynamics that involve 5-. 


0 = Mfj + Cfl + (K + qD Q )r\ + qC Q x 

= Mi) + [Cfl + (K + q 0 D Q )q + ? 0 C Q x] + 8-[D Q n + C Q x] 

= t) + [M 'Cli + M J (K + g 0 Dn)' r l + <7 0 M CqX] 

(42) 

+ 5^[M _1 D q ti + M _1 C Q x] 

= i\ + [M -1 Cfl + M" 1 (K + q 0 D Q )ri + <? 0 M -1 C q x] + S-z 
= f\ + [M -1 Cfi + M _1 (K + q 0 D Q )Tl + ^M _I C q x] + w 

The signals z and w are introduced into this formulation to associate the perturbation in dynamic pres- 
sure to the nominal dynamics in a feedback manner. The signal z can be generated as an output of the 
plant because z is a linear combination of states. 

z = M _1 DqT) + M _1 CqX (43) 
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The signal w is related to z by the dynamic pressure perturbation. 


w = 5- z (44) 

H 

The state-space aeroelastic model for nominal stability analysis in the p. framework is formulated 
using the state-update matrix. This matrix is determined by the dynamics at the nominal dynamic pres- 
sure, and the additional input and output signals that introduce perturbations to the dynamic pressure. 
That perturbation, 5- , is not an explicit parameter in the state-space model because 5- only affects the 
plant through a feedback relationship as determined by the signals z and w. Define the transfer function 
P(s) generated by state-space matrices such that z = P(s)w. 


11 


0 

I 

0 

0 

ii 


-M“‘(K + ? 0 D Q ) 

-M _1 c 

-q M _1 C q 

-I 

X 


b q 

0 

a q 

0 

Z 


[ M-'Dq 

0 

m-'c q 

0 


(45) 


Figure 4. 1 shows the feedback interconnection between the perturbation in dynamic pressure and the 
nominal plant model parameterized around that perturbation. This interconnection is an LFT, and the 
small gain condition of lemma 2.4.1 or the p condition of theorem 3.3.3 can be directly applied to analyze 
stability with respect to a variation in the flight condition q . 



Figure 4.1. Linear fractional transformation system for nominal stability analysis in the fi framework with 
parameterization around perturbation in dynamic pressure. 


Formulating the nominal aeroelastic dynamics in the |l framework immediately demonstrates the pro- 
cedure used in computing a flutter margin. Traditional flutter analysis algorithms such as the p-k method 
and the ft method as applied to figure 4.1 are searching for a value of dynamic pressure that results in a 
flutter instability. The nominal flutter margin question may be posed, which is answered by these methods. 

Question 4.2.1 (nominal flutter margin): What is the largest perturbation to dynamic pressure for 
which the nominal aeroelastic dynamics are stable? 

The dimension of the uncertainty block is the dimension of the signal z. The state-space equations for 
p (s) demonstrate this dimension is the number of modes in the system, n. The number of free variables in 
the ft upper-bound optimization, and consequently the computational cost of p,, is a function of the 
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uncertainty dimension. In this way, the number of aerodynamic states, n Q , does not directly affect the cost 
of the flutter estimation. The only cost increase caused by these additional states is computing the fre- 
quency response of the state-space matrix, which is generally much lower than the cost of computing p. 

Demonstrating the aeroelastic system formulated in the p framework in figure 4.1 is straightforward. 
The dynamic pressure parameterization is equivalent to the nominal state-space system given in t e pre- 
vious section. Simply compute the closed-loop transfer function with no uncertainty, 5- = 0, and the 

nominal system is recovered. 


Wind-tunnel and ground vibration testing can experimentally determine aerodynamic stiffness and 
damping matrices that are more accurate than the matrices approximated by a finite-element model. 
These matrices can be incorporated into a nominal state-space model in the p framework using the deri- 
vation in Appendix B. 


This procedure considers variations in dynamic pressure for an aeroelastic model at constant Mach 
number. The unsteady aerodynamics are highly nonlinear with variation in Mach number, ^and attempts to 
model Mach variations in an LFT may produce highly conservative flutter margins. ’ The p method 
presented herein is considering flutter margins in terms of dynamic pressure as measured along lines of 
constant Mach Flutter is a function of the two variables, dynamic pressure and Mach, so computing the 
dynamic pressure causing flutter for a dense set of discrete Mach values will generate an accurate portrait 

of the flutter margins. 


4.3 Robust Aeroelastic Model in the Structured Singular Value Framework 


A robust aeroelastic model in the p framework can be generated by associating uncertainty operators, 
A, with the nominal model and including the parameterization around a perturbation in dynamic pressure. 
These uncertainty operators can resemble any of the forms presented in section 3.2, including parametric 
uncertainty and additive and multiplicative representations of dynamic uncertainty. 

Choosing a reasonable uncertainty description is crucial for determining a valid robust flutter margin. 
This choice can arise logically from consideration of weaknesses in the modeling process, previous 
experience with aeroelastic analysis, and comparison with observed flight dynamics. Chapter 5 gives a 
noncomprehensive examination of several obvious uncertainty descriptions that may be associated with 

an aeroelastic model. 

The LFT is a valuable framework for formulating the robust aeroelastic model so that the model is 
suitable for p analysis. The various system blocks composed of the nominal state-space model with asso- 
ciated uncertainties and any additional subsystem blocks and their associated uncertainties can be 
expressed as a single model and uncertainty operator. 

Figure 4.2 shows the block structure used for p analysis of the uncertain aeroelastic system. The 
structured singular value is computed with respect to a single, block-diagonal, structured operator that 
contains the perturbation to dynamic pressure and the structured uncertainty operator along the diagonal. 
The perturbation to dynamic pressure is explicitly shown to distinguish 8- from the modeling uncer- 
tainty and emphasize 8- as a special operator used to describe a range of flight conditions. 
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Figure 4.2. Linear fractional transformation system for robust stability analysis in the p framework with 
parameterization around perturbation in dynamic pressure and structured uncertainty. 


The flutter margin computed for the uncertain system (fig. 4.2) is a more accurate margin than one 
computed with traditional methods such as p-k. These traditional methods address the nominal flutter 
problem in question 4.2.1; the robust flutter margin must consider the effect of the modeling uncertainty. 
The robust flutter margin actually finds the smallest perturbation to dynamic pressure for the entire set of 
plants formulated by the interconnection of the nominal dynamics and all elements A € A. Question 4.3.1 
poses how to compute these margins. 

Question 4.3.1 (robust flutter margin): What is the largest perturbation to dynamic pressure for 
which the nominal aeroelastic dynamics are robustly stable to the entire range of modeling errors as 
described by the norm-bounded A? 

This question can be answered by computing p for the system in figure 4.2. 

4.4 Computing a Flutter Margin with the Structured Singular Value 

The nominal flutter problem posed by question 4.2.1 and the robust flutter problem posed by 
question 4.3. 1 can be solved as a |ll computation. The value of p is a sufficient direct indication of the 
flutter margin for the nominal system; additional information regarding the norm bound on the uncer- 
tainty is required to derive the robust flutter margin. 

The nominal aeroelastic model is formulated for stability analysis in the |l framework in section 4.2 by 
introducing a perturbation, 8- , to dynamic pressure. A nominal flutter margin is computed to answer ques- 
tion 4.2.1 by considering the smallest value of this perturbation that destabilizes the model. The nominal 
flutter pressure can be directly calculated by computing p with respect to the perturbation operator 5- . 

The exact value of p can be analytically formulated to compute a nominal flutter pressure because a 
closed-form solution for p with respect to a single, real, scalar operator is known. This solution is the 
spectral radius of the frequency-varying transfer-function matrix. 

p(P) = max p(PO'G))) 

coeR (46) 

The spectral radius of P(yco) is a discontinuous function of frequency, so computational algorithms 
based on searches over a finite set of frequency points may not guarantee the correct computation of 
robustness values. A small amount of complex uncertainty can be added to the real uncertainty that 
allows a continuous p function to be analyzed but introduces unrealistic conservatism. An heuristic 
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robustness indicator can be substituted for p that considers stability over a frequency segment but is not 
34 35 

considered here. ’ 

A relatively simple approach can be used to compute p with respect to a single real parameter by 
considering the destabilizing value of the parameter. A search over the parameter space will result in 
computation of p and the desired flutter margin. Lemma 4.4. 1 presents the principle of this approach. 

Lemma 4.4.1: Given the plant, P, derived at nominal dynamic pressure, q 0 , with a perturbation to 
dynamic pressure, 5-, arranged in the feedback relationship of figure 4.2, define 5. 

8 = min { 8 ^ : F U (P, 8 -) is unstable} 

8->o q 

Q 

Then p(P) = l such that gutter = 4 0 + 5 is the nominal flutter pressure and 

T = 8 

nom 

represents the nominal flutter margin answering question 4.2.1. 


Proof: 

This result is immediately obvious using definition 3.3.1 for p with respect to a scalar uncertainty 
parameter 8- The p is the inverse of the smallest destabilizing perturbation, and 5 is computed as the 
smallest value of 8 5 that destabilizes P; thus p(P) = g . 

□ 

Lemma 4.4.1 indicates a computational strategy to compute a nominal flutter margin that does not 
require a search over a set of frequency points. The flutter pressure is found by increasing values of 8- 
until an eigenvalue of the state matrix of F U (P, 8- ) has a negative real part indicating F u (P, 8- ) is unsta- 
ble. Algorithm 4.4.2 demonstrates a bisection search implementation that efficiently computes upper and 
lower bounds on the minimum destabilizing 8- perturbation to within a desired level of accuracy. 

Algorithm 4.4.2 (nominal flutter margin): 

Given plant P at nominal dynamic pressure, q 0 , affected by perturbation 8- as in figure 4.1: 

Define scalars 8 upper > 8 lower >0 as bounds on the smallest destabilizing 8- . 

Define scalar e > 0 for accuracy. 

while U pp er ~ ^ lower '* ^ 

= 2 ^upper + ^loweri 

if F U (P, 8 -)hasan unstable pole), then § upper = 8^ 


else 


} lower 




_nom _ _ ^ 

^ flutter ~ + 0 upper 

y £ 

1 nom ° upper 
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Robust flutter margins that address question 4.3.1 cannot be computed using algorithm 4.4.2 because 
an additional search over the set of uncertainty operators A must be included. These margins must be 
computed using the augmented plant P, which includes feedback signals relating the perturbation to 
dynamic pressure and the uncertainty description as shown in figure 4.2. Define the block-structured set 
of operators, Ag_, which considers a particular perturbation 8- and set of operators A describing 
model uncertainty. 1 



\ ■ \ = 


0 

0 A 


-AeA, ||A|I ^ 1 


(47) 


A set of plant models, F (P,AgJ, exists for each value of 8- that defines the nominal plant at 
dynamic pressure q - q Q + 5- and variations to that plant caused by the set of uncertainties A. The 
robust flutter margin corresponds to the smallest perturbation 8- for which an unstable plant, F U (P,A g ), 
exists with Ag_ e A g_. If every member of the set of plants F(P,Ag ) is stable, then q = q + 8- is not 

q q * u v q 1 1 0 Cj 

a flutter pressure and P formulated at q is robustly stable to the uncertainty description A. 

The smallest destabilizing perturbation 8- corresponding to a robust flutter margin can be computed 
by a p computation. The (I framework analyzes robustness with respect to a single, structured operator, so 
the operator set used to compute a robust flutter margin must contain the set of uncertainty operators A 
along with a range of dynamic pressure perturbations. The p will compute the robustness of P with 
respect to this operator set to find the smallest destabilizing perturbation to dynamic pressure and the 
smallest destabilizing uncertainty operator. Define the A that contains Ag_ sets for a norm-bounded set 
of 8- operators. * 


A = 


f T - 

8- o" 

! A : A = 

_o aJ 


, Ae 


A, HAIL <1, ||8-| 



(48) 


Imposing the norm bound for 8- operators as ||8-|^< 1 may seem overly restrictive because the 
units of 8- are the same as the units of q in the model. This condition implies the A considers the range 
of flight conditions q =q o ± 1 lbf/ft 2 for plants formulated by dynamic pressure in units of lbf/ft 2 . Such 
a small range of flight conditions is not useful for stability analysis unless q o is extremely close to the 
flutter pressure. This limitation is avoided by introducing a weighting function, W- , to the computation 
of q . q 


q = q 0 + (49) 

A W- >1 allows a large range of flight conditions to be considered despite the unity norm-bound 
constraint on 8- . This weighting is incorporated into the stability analysis by scaling the feedback signals 
between the 8- operator and the plant P to form the scaled plant, P . 
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p = p 


W-0 
0 I 


(50) 


A robust flutter margin is computed by analyzing p(P) with respect to the_ A. The r^ust flutter 
nressure is determined by iterating over scalings W- until the smallest pressure q ^o + q 
KSi, Stable to the set of uncertainties A. Theorem 4.4.3 formally demonstrates 

this concept. 

Theorem 4.4.3: Given the plant P derived at nominal dynamic pressure q Q with a perturbation to 
dynamic pressure b- a and set of uncertainty operators A norm-bounded by one arranged in the feedback 
relationship of figure 4.2, define the plant P with real diagonal matrix W- scaling the feedback signals 

relating 5- and P. 


P = P 


W-0 
0 I 


(51) 


Then q^^r - 4 0 + is the robust flutter pressure if and only i/p.(P) - 1. Also, 


r - = W- 

rob q 


(52) 


represents the least conservative robust flutter margin answering question 4.3.1. 


Proof 


<= (necessary) 

The condition p(P) = 1 implies that the smallest destabilizing perturbation to P is described by some 
A e A with || a|| ^ = 1 , so no destabilizing A e A exists and the smallest positive destabilizing perturba- 
tion to dynamic pressure is at least 8- = 1 , which corresponds to dynamic pressure q = q 0 +W-b-- 
q + w_ Thus, P is guaranteed to be robustly stable to the uncertainty set A for any perturbation to 

dynamic pressure less than W-, so W- is a robust flutter margin. 

=> (sufficient) 

Assume p(P) > 1. Define real, scalar a < 1 such that p(P) = 5 , which implies, from theorem 3.3.2, 
that P is robustly stable to all uncertainties A 6 A with 0L<a< 1 . Thus, P is not guaranteed to be 
stable over the entire range of modeling uncertainty defined by the unity norm-bounded set A, so this 
perturbation is not a valid robust flutter margin and does not answer question 4.3.1. 

Assume p(P) < 1. Define real, scalar a > 1 such that p(P) = 5 , which implies, from theorem 3.3.2, 
that P is robustly stable to all uncertainties Ae A with |aL<cc. Thus, P is robustly stable to an 
uncertainty description larger than that defined by the unity norm-bounded set A, so this condition 
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defines a valid flutter margin but is not the least conservative robust flutter margin and does not answer 
question 4.3.1. 

□ 


The only difference between models P and P results from the scaling W- , which scales the feedback 
signals between P and 5-. No external scaling matrix is allowed to affect the feedback signals between 

H 

P and A because A is defined with a unity norm bound. Computing p of the plant P with an additional 
scaling on the lower-loop signals would consider a scaled set of operators A that does not accurately rep- 
resent the uncertainty description. Therefore, P only scales the 5- feedback signals. 

Theorem 4.4.3 can be modified to compute a nominal flutter margin by changing the identity matrix 
in the scaling used to compute P to a zero matrix. This modification eliminates the feedback interconnec- 
tion between the model and the uncertainty description A, so p considers only the nominal dynamics and 
computes the smallest destabilizing perturbation to dynamic pressure and T rob = T . 

Including the uncertainty description A ensures the robust flutter margin will be no greater than the 
nominal flutter margin. The robust flutter margin considers the model used to compute the nominal flutter 
margin that corresponds to the uncertainty operator A = 0 e A and the models that correspond to the 
remaining operators A e A . The conservatism in the robust flutter margin makes intuitive sense because 
the nominal flutter margin is the worst-case stability boundary for a single model and the robust flutter 
margin is the worst-case stability boundary for a family of models. 


T 


rob 


<r 

nom 


(53) 


The proof demonstrating the necessary and sufficient condition p(P) = 1 also makes intuitive sense 
because the equality sign is needed to ensure the flutter margin is valid without being overly conserva- 
tive. If p(P) < 1, then no AeA causes an instability and the flutter margin is too conservative. If 
|i(P) > 1, then the system is not robust to all modeling errors AeA and the flutter margin is not a valid 
robust flutter margin. 


Theorem 4.4.3_demonstrates a robust flutter margin can be computed by determining a scaling matrix 
W- for which |X(P) = 1. Algorithm 4.4.4 implements an iterative approach to compute a scaling matrix 
for which p(P ) e (1 ± e) for some desired level of accuracy e. 

Algorithm 4.4.4 (robust flutter margin): 


Given plant P at nominal dynamic pressure q affected by unity norm-bounded 8 - and A as in 
figure 4.2: q 

Define initial weighting W-. 

H 

Define scalar £ > 0 for accuracy. 


P = P 


0 


0 

ij 
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while (fl(P) > 1 + e) OR (fl(P) < 1 - e) { 

W, 

W- = 


p(P) 
P = P 


W-0 


0 I 


} 

iflLr = % + W < 

r TO „ = w, 


The dynamic pressure defining the flight condition for the nominal plant dynamics must be chosen 
carefully to ensure the robust flutter margin computed with algorithm 4.4.4 is valid. The nature of fl an 
the upper bound is such that all norm-bounded operators centered around the origin are assumed to be 
valid perturbations to the plant dynamics, so positive and negative perturbations to dynamic_pressure are 
considered by the robust stability analysis. Thus, the robust flutter margin computed by p(P) - 1 could 
correspond to either perturbation 8^ = 1 or 8^ — —1. 

The actual dynamic pressure at the flutter instability is qj^ter = lf h = A large weightmg 

of W- > q indicates the flutter occurs at a negative dynamic pressure that may be unrealistic for classical 
flutter analysis. The value of the nominal dynamic pressure q 0 can slide along the real axis to a large value 
without loss of generality in the p analysis because this parameter linearly affects the nominal dynamics. 
A simple approach to ensure that the robust flutter pressure is a positive dynamic pressure is demonstrated 
in algorithm 4.4.5, which iterates over increasing values of q Q until the scaling associated with the robust 

flutter margin satisfies W - < q 0 . 

Algorithm 4.4.5 (robust flutter margin with q 0 iteration): 

Given parameters as in algorithm 4.4.4: 

Given initial value of nominal dynamic pressure q Q \ 

valid_margin = FALSE 

while (validjnargin == FALSE) { 

compute plant P at nominal dynamic pressure q 0 

compute q fiu tter and associated W - from algorithm 4.4.4 

if(v- q >q 0 ^ then 4 o = LI 

else validjnargin = TRUE 

} 

The exclusion of low and negative dynamic pressures for stability analysis may not be desirable for 
all applications related to aeroelasticity. An example of such an application is the analysis of aeroser- 
voelastic dynamics for a high-angle-of-attack aircraft that concerns instabilities at low dynamic 
pressures. 36 The procedure for these types of analysis chooses a low value of q 0 and finds the scaling 
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W- corresponding to 5- = -1, which computes the low dynamic pressure instability. Algorithm 4.4.5 
can be modified for these applications. 

The flutter margin computation must allow for an arbitrary structure of operators A, so an upper 
bound such as the function derived in appendix A must be used. Algorithms 4.4.4 and 4.4.5 can be 
adapted by replacing the |i calculation with a |i upper-bound calculation to compute flutter margins. A 
search over frequency points is required when using the upper bound (appendix A), so the accuracy of the 
robust flutter margin requires a dense grid of frequencies associated with the natural frequencies of the 
worst-case dynamics to be considered. 

A simple approach can be implemented if the natural frequency of the unstable dynamics at the nom- 
inal flutter pressure can be assumed to be similar to the natural frequency of the unstable dynamics at the 
robust flutter pressure. This assumption can be justified if the uncertainty does not change the critical 
flutter mode between the nominal and robust flutter pressures, which is often true for systems that have 
relatively small levels of uncertainty and clear separation between critical and subcritical modal frequen- 
cies. Algorithm 4.4.6 presents this approach, which first computes the frequency of the nominal flutter 
mode and then computes a robust flutter margin from the p upper bound evaluated near that frequency. 

Algorithm 4.4.6 (robust flutter margin with reduced frequency grid): 

Given the system in figure 4.2: 

Compute frequency to associated with nominal flutter dynamics using algorithm 4.4.2. 

Define dense frequency grid Q centered around (0. 

Compute robust flutter margin from p upper bound evaluated at Q using algorithm 4.4.5. 

A large and dense frequency grid increases the confidence that computed robustness measure is actu- 
ally an upper bound for p. Algorithm 4.4.6 must be used with caution because the assumptions behind its 
use may not be satisfied. Computing a robust flutter margin with a dense frequency grid for a particular 
aircraft at several Mach numbers and then comparing the frequencies of the flutter dynamics to those of 
the nominal flutter dynamics is recommended. If these frequencies are similar, then algorithm 4.4.6 can 
be considered for further analysis at different Mach numbers. 

4.5 Properties of the Structured Singular Value as a Flutter Margin 

The flutter computation method described herein uses p as the worst-case flutter parameter. The 
structured singular value is a much more informative flutter margin than traditional parameters such as 
pole location and modal damping, so p presents several advantages as a flutter parameter. 

The conservatism introduced by considering the worst-case uncertainty perturbation can be interpreted 
as a measure of sensitivity. Robust p values that are significantly different than the nominal flutter margins 
indicate the plant is highly sensitive to modeling errors and changes in flight condition. A small perturba- 
tion to the system can drastically alter the flutter stability properties. Conversely, similarity between the 
robust and nominal flutter margins indicates the aircraft is not highly sensitive to small perturbations. 

Robustness analysis determines not only the norm of the smallest destabilizing perturbation but also 
the direction. This information relates exact perturbations for which the system is particularly sensitive. 
Thus, p can indicate the worst-case flutter mechanism, which may naturally extend to indicate active and 
passive control strategies for flutter suppression. 
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Damping is only truly informative at the point of instability because stable damping at a given flight 
condition does not necessarily indicate an increase in dynamic pressure will be a stable flight condition. 
The structured singular value computes the smallest destabilizing perturbation, which indicates t e 
nearest flight conditions that will cause a flutter instability. In this respect, p is a stability predictor and 
damping is merely a stability indicator. 

These characteristics of p make the worst-case flutter algorithm especially valuable for flight test 
programs. Aeroelastic flight data can be measured at a stable flight condition and used to evaluate 
uncertainty operators. Unlike damping estimation, the p method does not require the aircraft to approach 
instability for accurate prediction. The p can be computed to update the stability margins with respect to 
the new uncertainty levels. The worst-case stability margin then indicates what flight conditions can be 

safely considered. 

Safe and efficient expansion of the flight envelope can be performed using an on-line implementation 
of the worst-case stability estimation algorithm. Computing p does not introduce an excessive computa- 
tional burden because each F/A-18 flutter margin presented herein was derived in less than 2 min using 
standard off-the-shelf hardware and software packages. The predictive nature of p and the computational 
efficiency allow a flutterometer tool to be developed that tracks the flutter margin during a flight test. 
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CHAPTER 5 

UNCERTAINTY DESCRIPTIONS IN AEROELASTIC MODELS 


5.1 Parametric Uncertainty in Structural Models 

Parametric uncertainty denotes operators that describe errors and unmodeled variations in specific 
elements and coefficients in dynamic system equations. Recall the generalized aeroelastic equation of 
motion for state vector r\ e R n , 21 


Mrj + Cf] + Kr| +^Q(s)ri = 0 (54) 

Robust flutter margins computed with the |i method are strongly affected by the choice of uncertainty 
descriptions associated with these dynamics, so this uncertainty must be a reasonable indicator of poten- 
tial modeling errors. Parametric uncertainty can be directly associated with the structural matrices to indi- 
cate specific errors in the finite-element model. 

Define an operator, A K € R n x n , that describes additive uncertainty of a nominal stiffness matrix 
K Q . Associate a weighting matrix, W R e R n n , with this uncertainty such that a stability analysis 
should consider a range of stiffness matrices described by all A K with ||A k | < 1 . 

K = K o + W k A k (55) 

Parametric uncertainty can also be associated with structural elements in a multiplicative relationship. 
Define an operator, A c e R n , that describes multiplicative uncertainty of the nominal damping 
matrix C Q . A weighting W c e R "is again associated with the uncertainty such that the anticipated 
range of damping matrices for robust stability analysis is described by all A c with |A C || < 1 . 

C = C 0 (I + W c A c ) (56) 

The choice of additive uncertainty for A K and multiplicative uncertainty for A c does not reflect any 
generalized assumptions regarding the proper way to model errors in stiffness and damping; rather, each 
type is included to demonstrate the different mathematical derivations. Additive and multiplicative 
operators are common types of uncertainty models, so demonstrating how these types of uncertainty are 
associated with a structural model is instructive. The actual choice of A^ and A^. is problem-dependent 
and can vary with different aircraft and different finite-element modeling procedures. 

The uncertainty operators are described in this section as A R , A c e R n x n w i t h real elements 
because the operators describe perturbations to the generalized stiffness and damping matrices that are 
usually real. These operators are often additionally constrained to be diagonal operators with n indepen- 
dent parameters because the generalized stiffness and damping are often computed as real diagonal 
matrices. The real and diagonal nature of these uncertainties is not required for p analysis, so full-block 
complex uncertainties can be used if they better describe the nature of the modeling errors. 

Also, the weighting functions W K , W c € R nX n are presented as constant and real matrices because 
the functions are associated with constant and real stiffness and damping matrices. These constraints on 
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the weightings can be relaxed if the nature of the uncertainty is best described by complex and frequency- 
varying weighting functions. 

Substitute the uncertain K and C into the equation of motion, including the state-space representation 
of the unsteady aerodynamic forces Q(s) presented in section 4.1. Introduce a perturbation, S-, to 
dynamic pressure, and separate the nominal dynamics from the unknown terms. 

0 = Mf| + Cfi(K + <?D q )t) + <?C Q x 

= Mri + c 0 (I + W c A c )li + (K 0 + W k A k + ( q Q + 8^)D q )ti + ( q a + 5^)C Q x 
= ri + [M -1 C 0 fi + M 1 (K 0 + q a Dq )r| + q Q M CqX] 

+ 5-[M' 1 D Q r| + M“ 1 C Q x] + A K [M' 1 W K Ti] + A c [M ‘c^Weti] 

= fi + [M _1 C 0 f] + M" 1 (K 0 + q Q Dq )T| + q 0 M _1 C Q x] + 8- z- + A K z R + A c z c 
= ri + [M' 1 C 0 f] + M' 1 (K 0 + q 0 D q )ti + q 0 M" 1 C Q x] + w- + w K + w c 

The signals z- and w- are introduced in section 4.2 to relate the perturbation in dynamic pressure to 
the nominal plant through a feedback relationship. The feedback operation w- = b- z- is used where z- 
is a linear combination of the states of the plant. 


z- = M *D q T| +M ‘c q x (58) 

Additional signals are introduced to the plant formulation, where z R and w R are associated with the 
uncertainty in the stiffness matrix, and z c and w c are associated with the uncertainty in the damping 
matrix. The outputs of the plant z R and z c are formulated as linear combinations of the states. 


z R = M l W R r| 
z c = M-‘w c C 0 f| 


(59) 


The feedback mechanism to describe the modeling uncertainty uses a relationship between these out- 
put signals and the w R and w q input signals. 


w R = A r z r 
W C = A C Z C 


(60) 


The state-space plant matrix can be formulated using these additional input and output signals. 
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Figure 5. 1 shows how the uncertainty operators and perturbation to dynamic pressure affect this plant 
formulation in a feedback manner. 



Figure 5.1. Linear fractional transformation system for robust stability analysis in the fi framework with 
parameterization around perturbation in dynamic pressure and uncertainty in structural stiffness and 
damping matrices. 


The formulation does not directly allow uncertainty in mass to be described by a feedback operator. 
The LFT and |i frameworks require uncertainty operators to affect the nominal dynamics in a linear man- 
ner, and this requirement precludes introducing mass uncertainty. The inverse of the mass matrix scales 
most terms in the state matrix of P(s), including terms involving q. Associating a mass uncertainty oper- 
ator A m with the mass matrix scaling q would introduce terms of A M 5- , which is a nonlinear function 
of uncertainty operators and cannot be directly considered by the p method. 

5.2 Parametric Uncertainty in Aerodynamic Models 

The unsteady aerodynamic forces Q(s) e C n can be represented as a state-space model with 
nq states. 

Q(s) = D Q + C Q (si - A Q f ' b q (62) 
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Parametric uncertainty can be associated with the matrix elements of this state-space representation to 
describe errors. These errors can result from several sources in the modeling procedure including 
computational fluid dynamic algorithms that determine the frequency-varying forces and the system 
identification algorithms that represent the computational values as a slale-space system. 

n n xn, 


Define an operator A Aq e R 


1qX ° Q to describe uncertainty in the state matrix of Q(s). This operator 
directly affects a nominal^ and describes errors and variations in the poles of tte state-space represen- 
tation of the unsteady aerodynamic forces. Include a weighting function W Aq e R such that the range 
of state matrices to be considered by robust stability analysis is described by all A Aq with || A A( II < 1 • 


A Q = a q 0 ' r VY A ~^ A 


+ W A A a 
Q A Q 


(63) 


Define also an operator A B e r" 0 *" to describe multiplicative uncertainty in a nominal B Qq 
matrix of Q(s). A weighting function W Bq eR Q is associated with this uncertainty such that the 
range of possible B Q matrices is described by all A Bq with |A B ^ 1. 

Bq = B Qo (I + W Bq A B q ) m 


The choice of additive and multiplicative operators is made for reasons similar to those presented in 
section 5.1. One of each type of uncertainty is included to demonstrate the derivation procedures o ow 
each uncertainty is associated with the nominal aeroelastic dynamics in a feedback relationship. The 
actual choice of which type of uncertainty is most suitable to describe errors in A Q and B Q is 

problem-dependent. 

Also, defining the uncertainty operators as real and weightings that are real and constant is not a 
requirement. This section presents the problem formulation with these definitions because associating 
these types of uncertainties and weightings with the constant real A Q and B Q matrices makes intuitive 
sense. Certainly p can also compute robust flutter margins with respect to complex frequency-varying 
weighted uncertainties. 

The nominal aeroelastic model in section 4.2 defined the vector x as the states associated with the 
state-space model Q(s). The matrices A Q and B Q directly affect the aeroelastic dynamics only through 
the state derivative equation for x and do not appear in the other state derivative equations. The aeroelas- 
tic plant can be formulated in the p framework by substituting the uncertain values of A Q and B Q into 
this state derivative equation without considering changes in the derivative equations for the remaining 

states rj and p . 


x = 


AqX + BqTI 


= (V + W A Q % )X + B Q„ (I + W B Q A B Q )n 
- a Q„ x + b Q o 11 + A a 0 w a 0 x + a b q w b q b q o 1 
= a q„ x + b Q 1 , 11 + A a„ z a.. + a b, z b 


Q Q 


Q U Q 


= A Q o x + B Q o n + w A Q + w BQ 


( 65 ) 
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Several signals are introduced to this equation, where z A and w Aq are ass o c i a ted with the uncer- 
tainty in Aq and z B , and w B is associated with the uncertainty in Bq. The signals z Aq and z B ^ are 
output from the plant matrix as linear combinations of the states. 


Z A = w a x 
a q a q 

Zd — Wn Br, T| 

B 0 B Q 


( 66 ) 


The feedback mechanism to describe the modeling uncertainty uses a relationship between these 
output signals and the w K and w c input signals. 


w A = A a z a 
a q a q a q 

w R = A r z r 
b q b q b q 


(67) 


The state-space plant matrix can be formulated using these additional input and output signals. 
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Figure 5.2 shows how the uncertainty operators and perturbation to dynamic pressure affect this plant 
formulation in a feedback manner. _ 



Figure 5.2. Linear fractional transformation system for robust stability analysis in the |i framework with 
parameterization around perturbation in dynamic pressure and uncertainty in Aq and Bq matrices of the 
state-space Q(s) model. 
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Associating a A B uncertainty operator with the B Q matrix may not seem immediately useful 
because considering errors in the poles determined by the A Q matrix is often intuited. This A Bq uncer- 
tainty can be essential to accurately describe the modeling errors because errors in terms common to both 
Aq and Bq may exist. Such a situation arises for certain modeling representations of aerodynamic lags. 
Consider a simplified Roger’s form Q(s) matrix that uses two Pade approximations to represent lag terms. 



-P, 0 

-Pi 

Q(s) = — V + -^o" = 

S + Pj S + P 2 

0 -p 2 

-P 2 

1 1 

2 


(69) 


The poles of this system are determined entirely by the Aq matrix, so uncertainty in the poles can 
be entirely described by a A^ operator associated with the Aq matrix. A similar uncertainty A B ^ 
should also be associated with the B Q matrix in this case because the Pj and P 2 terms appear in both Aq 
and Bq. Allowing variation in Aq but not in Bq will introduce unwanted zeros to the system, so the 
proper way to model pole uncertainty for this formulation is to include both A^^ and A B ^ °P era ^ ors - 

The Padd approximation appears often in aeroelastic models, so demonstrating the LFT formulation 
of Q(s), which includes the uncertainties in poles p j and P 2 in a feedback relationship, may be useful. 
The uncertainties in Q(s) can be developed distinctly from the structural uncertainties because LFT oper- 
ations allow the structural and aerodynamic models to be combined into a single plant model with a struc- 
tured uncertainty description. 10 


Define real scalar operators A B and A B to describe uncertainty associated with nominal values of 
the poles pj and P 2 . Real scalar weightings Wp and Wp normalize the uncertainty such that the 
range of pole°s to be considered by robust stability analysis is described by all A B ^ with A B |^ < 1 and 


all Ap^ with 


J 2 ~ 


< 1. 



(70) 


Define states and x 2 of Q(s), and consider an input signal Uq that generates the output signal yQ 
by the relationship yQ = Quq. Substitute the uncertain Pj into the state equation of x } . 


Xj - “PjXj - P i U Q 

= -(P u + W p A P| )x,-CP 1o + w p, a p, )u q 

= -Pi/i-Pi^q-V^i+w^uq) (71) 

= -Pi„ x i-Pi 0 u q-Vp, 

' -P| 0 X I -Pl 0 u Q - W P, 
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Perform a similar derivation for the state equation of x 2 . 


x 2 - P 2 x 2“P2 u Q 

= -(P 2o + Wp 2 Ap 2 )x 2 -(P 2o + Wp 2 Ap 2 )u Q 

= -P 2 ( x 2 - P 2 o U Q - Ap 2 (Wp 2 X 2 + Wp 2 U Q ) (72) 

= -P 2 q X 2 - P 2 ( U Q - Ap^ 

= — P 2q x 2 — ^2 0 U Q — w p 2 

The signals Zp ( and w p are introduced as plant output and input signals to related the uncertainty 
^Bj * n a feedback manner. The signals Zp^ and Wp^ are similarly introduced to relate the uncertainty 
A b 2 in a feedback manner. The state-space matrix can be formulated to describe the nominal Q (s) 
with these additional input and output signals. 
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Figure 5.3 shows the block diagram for Q(s) and the uncertainties. 
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Figure 5.3. Linear fractional transformation system describing Pade approximation to represent unsteady 
aerodynamic force matrix in the jn framework with uncertainty in lag terms. 


5.3 Dynamic Uncertainty 

Dynamic uncertainty operators are often associated with aeroelastic models to account for modeling 
errors that are not efficiently described by parametric uncertainty. Unmodeled dynamics and inaccurate 
mode shapes are examples of modeling errors that can be described with less conservatism by dynamic 
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uncertainties than with parametric uncertainties. These dynamic uncertainties are typically complex 
order to represent errors in both magnitude and phase of signals. 

Consider a system P having two modes with natural frequencies at 4 rad/sec and 30 rad/sec. 


in 


P = 


16 1 

900 | 

,S 2 + 0.4s + 16y 

^s 2 + 0 .6 + 900 J 


(74) 


Define a nominal model P„, which will be used for stability analysis of the system but does not 
model the high-frequency mode of P. 


P 0 = 


16 


\ 


U 2 + 0.4s + 16 


(75) 


The large difference in natural frequency between the high- and low-frequency modes of P precludes 
par^etn” uncertainty associated with the low-frequency mode of P D from be.ng a reasonab 
description of the modeling errors in P„. The magnitude of any parametric uncertainty associated with 
die low-frequency mode would need to be extremely large to account for the unmodeled high-frequency 
dynamics, so the stability analysis would be relatively meaningless because the large uncertainty wou 
imply the plant is not accurate at any frequencies. 

. \ ez C' fan he used to describe the high-frequency modeling 

varying levels of modeling errors. 


W = 100 


s + 0.5 


_ (76) 

s + 500 

The set of plants P used for robust stability analysis is formulated to account for the range of dynam- 
ics as described by the norm-bounded multiplicative uncertainty A. 


I JlClUlli VJ **•' 

l multiplicative uncertainty A. 

p = {P Q (I + WA) : ||AL<1} 

Figure 5.4 shows the block diagram for robust stability analysis of P, 


(77) 



Figure 5.4. Family of plants P = P 0 (I + AW) with input multiplicative uncertainty. 


42 





Figure 5.5 shows the magnitude of the transfer function from input to output for P and P , and the 
maximum magnitude of |P Q (I + WA)| , which is an upper bound for the output of P at each frequency. The 
multiplicative uncertainty is able to bound the modeling error at the high-frequency mode without intro- 
ducing excessive conservatism from large uncertainty associated with the low-frequency mode. 



Figure 5.5. Transfer functions for example system with multiplicative uncertainty. 


Dynamic additive operators may also be required in the uncertainty description to account for errors 
that are not efficiently described by either multiplicative and parametric uncertainties. Modeling errors 
associated with a zero of the system dynamics are an example of an error that is best described by addi- 
tive uncertainty. Multiplicative operators are not useful in this case because P Q (j to) = 0 at the 
frequency co associated with the zero of nominal model and, correspondingly, every member of the set of 
plants P Q 0' G))(I + WA) = 0 at this frequency. Additive uncertainty allows the system output for some 
member of the set of plants to be nonzero even at frequencies of the zeros of the nominal plant. 

Consider a plant P with several poles and zeros. 


P 


16 

64 


fs 2 + 0.48s + 64^1 

900 1 

+ 0.48s + 16y 

+ 0.6 s + 900 j 


(78) 


Assume the nominal plant P Q of this plant is similar to P and has the correct number of poles and 
zeros, but the coefficients of the system equations are incorrect. 


16 

fs 2 + 0.36s + 36 

900 ) 

36 

vOs 2 + 0.4s + 16 j 

k$ 2 + 0.6s + 900 j 


(79) 


43 



Define an additive uncertainty operator A e C that is normalized by a complex, scalar function W to 
reflect the frequency-varying levels of the modeling errors. 


W = 0.1 


s + 50 
s + 5 


(80) 


The set of plants 7> used for robust stability analysis is formulated to account for the range of dynam- 
ics as described by the norm-bounded additive uncertainty A. 


V= {P Q + WA : ||A|| M < 1} 

Figure 5.6 shows the block diagram for robust stability analysis of V. 


(81) 



Figure 5.6. Family of plants P=P 0 +WA with additive uncertainty. 

Figure 5.7 shows the magnitude of the transfer function from input to output for P and P Q , and the 
maximum magnitude of IP + WAl at each frequency. The additive uncertainty bounds the modeling 
error at each frequency, including the frequencies near the zero of the nominal plant, because the output 
of P is bounded above by the maximum magnitude of the members of the set V. 


Magnitude 



Figure 5.7. Transfer functions for example system with additive uncertainty. 
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These multiplicative and additive uncertainties are particularly important when comparing an analyti- 
cal transfer function with experimental transfer functions derived from flight data. Analytical models are 
often computed for a low range of frequencies because the high frequencies add complexity to the model 
but do not always affect the stability margins of the aircraft. The experimental data may indicate a high- 
frequency mode that is not included in the analytical model, so a frequency- weighted dynamic multipli- 
cative uncertainty can be associated with the model. 

The issue of mode shape uncertainty is often encountered when comparing low-frequency-predicted 
dynamical responses with flight data because sensor measurements are directly affected by the mode 
shapes. Both multiplicative and additive uncertainties may be required to accurately model mode shape 
errors and account for inaccurate response levels (which may be higher than predicted at some frequen- 
cies but lower than predicted at others) and inaccurate frequencies associated with poles and zeros. 

5.4 Uncertainty Associated with Nonlinearities 

The |X framework described uses linear operators to represent dynamical models and associated 
uncertainties but does not admit nonlinear operators. The |Li framework is useful for analyzing aircraft 
stability despite the constraints of linearity because physical systems, which are always nonlinear, can 
often be approximated by linear models to a high degree of accuracy. A classic example of this situation 
notes the linearized dynamics are often an acceptable representation of an aircraft operating near trim 
flight conditions, so linear models work well in practice for control synthesis and stability analysis. 39 

Nonlinear dynamics cannot always be accurately described by a linearized dynamics model, so the 
stability analysis should consider the effects of these nonlineanties/* 9 The (j. framework can associate 
linear uncertainty operators with linear models to describe the errors resulting from some types of 
unmodeled nonlinear dynamics. The uncertainty does not actually represent the nonlinearity; rather, the 
uncertainty allows the linear system responses to vary with sufficient magnitude to bound a range of non- 
linear system responses. 

Separating the nonlinear dynamics that cannot be linearized from the nonlinear dynamics that can be 
accurately represented by linear models is useful. Separate uncertainty descriptions can be formulated for 
each dynamical block, and the resulting operators can be combined using the LFT framework to formu- 
late a single, linear, plant model with a structured uncertainty description. Figure 5.8 shows an example 
LFT system representation for a nominal plant P o and associated uncertainty A, and a nominal linear 
model N q and associated uncertainty representing a system element with nonlinear dynamics. 



970987 


Figure 5.8. Linear fractional transformation system with nominal models and associated uncertainty 
operators. J 
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The system shown in figure 5.8 is commonly used to describe the coupling between the aeroelastic 
. y , t affecting an aircraft through control surfaces. Actuators can display many types 
of nmdine^tjehavtf^ 8 ^d e should I bTc»^idered fn the aeroelastic stability analysis, because pilot and 
autopilot commands that maintain trim during flight ensure the control surfaces are , con ^°^ y 
movhic 1 The errors in linear models resulting from unmodeled actuator dynamics such as nonlinear 
stiffness pa^etTor hys.eres.s functions can somettmes be described by a linear uncerratnty operator. 

Consider the response y of a nonlinear system N that models a system that has a 
corresponding to a hardening spring that is valid for the bounded mpu s.gnal u e [-10. 10]. Such 
systein can represent an element of an actuator model or a nonlinear structural mo . 

(82) 


y = Nu = 2u + 0.02u + 0.0082u 


Define a linear nominal model N„ such that y = N 0 u approximates the response of N. 


N 0 = 2 


(83) 


Associate an additive weighting operator A N with N„ such that stability analysis constders the set of 
plants Tt. 

^={N 0 +A n : ||A n L<1} (84) 

Figure 5 9 shows that the maximum and minimum magnitudes of the responses of the set 71 are able 
to hom^Ae response of the nonlinear system N for the range u e [-10, 10]. These bounds are overly 
conservative throughout this operating range, but they achieve the desired goal of describing errors in t e 
linear system response resulting from the unmodeled nonlinearity. 



Figure 5.9. System responses for hardening spring example. 
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A similar procedure can be used to describe error caused by an unmodeled nonlinear softening spring. 
Consider the response y of a system represented by N that is valid for an input signal u e [-10, 10]. 

y = Nu = 2u + 0.005u 2 -0.005u 3 (85) 

Define a linear nominal model N Q such that y = N o u approximates the response of N, 

N 0 = 2 (86) 

Associate an additive weighting operator A N with N 0 such that stability analysis considers the set of 
plants “H. 

^ = {N o +A N : ||A n (_ < 1 } (87) 

Figure 5. 10 shows that the maximum and minimum magnitudes of the responses of the set H are able 
to bound the response of the nonlinear system N for the range u g [—10, 10], These bounds are also overly 
conservative throughout this operating range, but they achieve the desired goal of describing errors in the 
linear system response resulting from the unmodeled nonlinearity. 



Figure 5.10. System responses for softening spring example. 


Another nonlinearity that commonly affects aeroelastic systems is hysteresis. The response of a 
hysteretic system depends on the trend of the input signal such that an increasing input signal generates a 
certain response, but a decreasing input signal generates a different response. Such hysteresis dynamics 
are difficult to express as a simple mathematical formula, so for illustrative purposes, assume N is a 
nonlinear hysteresis function whose response y depends on the trend of the input signal. Define a linear 
nominal model N 0 whose response y= N Q u approximates the response of the hysteretic N. 
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( 88 ) 


N 0 = 2 


Associate an 
plants Tt. 


additive weighting operator A N with N„ such that stability analysis considers the set of 

*={N c +A n : M_S1> (89) 

Figure 5 1 1 shows that the maximum and minimum magnitudes of the responses of the set * are able 

hysteresis nonlinearity. 



Figure 5.11. System responses for hysteresis example. 


Explicitly constraining the operating region of the input signal u € [-10 10] can be important to 
devdopiiig rea^nable uncertainties to describe errors resulting from unmodeled nonlmeari Ues^e 
errors fn the linear model can grow excessively large when considering a large range of inputs, so the 
uncertainty magnitude would need to also grow excessively large. The conservatism resulting from such 
a large uncertainty description may be unacceptable and require the input range to be constraine o more 

reasonable limits. 

The uncertainty description, even for a constrained operating region, will usually be overly 
tive^vhen^escribhig modeling errors for certain parts of the operating region. The uncertainty is able to 
bound the errors in § figures 5.9, 5.10, and 5.11, but the maximum and minimum bound are clearly not 
optimal Some amount of conservatism is expected when describing errors resulting from unmodeled 
nonlinearities because a linear model, whether a single plant or a set of plants, will usually not be an 
accurate representation of a nonlinear system that cannot be linearized. 
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Also, several types of unmodeled nonlinearities exist that are frequently encountered in aircraft 
systems but are not easily described by linear uncertainty operators associated with the linear models. 
Examples of these types of nonlinearities include free play, dead-band responses, friction, and rate limit- 
ing of actuators. 

The robust stability margins computed from |i with respect to the linear uncertainty operators describ- 
ing unmodeled nonlinear dynamics will always be somewhat suspect, because this approach cannot 
consider stability properties unique to nonlinear systems such as bifurcation points and limit-cycle behav- 
iors. 42 ’ 43 This approach limits the usefulness of the |l method to systems for which the nonlinearities have 
small effects on the response and do not introduce nonlinear instabilities to the critical flutter mechanism. 

5.5 Uncertainty Associated With Flight Data 

A theoretical model with an associated uncertainty description can be an accurate representation of 
the aeroelastic dynamics of an aircraft, but responses from that model may not identically match flight 
data. Additional uncertainties can be associated with the model to describe errors that are observed 
between the predicted responses and the measured responses from a commanded excitation to the air- 
craft. These uncertainties do not necessarily indicate errors in the model; rather, these uncertainties indi- 
cate errors in the process used to generate aeroelastic responses and measure flight data. 

One source of error is an incorrect assumption of excitation force used to generate the predicted and 
measured responses. The measured excitation force associated with the flight data may not correctly 
account for poor hardware performance and nonuniform spectral distribution of the force. Also, inexactly 
phased excitation between multiple force mechanisms can excite modes that are not anticipated by a 
theoretical analysis. A frequency-varying dynamic uncertainty can be associated with the force input of 
the analytical model to describe errors in the excitation. 

The phenomenon of nonrepeatibility can cause discrepancies between predicted and measured 
responses from multiple occurrences of excitation signals. Nonrepeatibility affects flight data by intro- 
ducing slight variations in responses, even for data recorded at identical flight conditions with identical 
excitation signals. This unexplained behavior may result from some unmodeled nonlinear dynamic or 
inexact excitation that is not correctly measured. External disturbances such as wind gusts or turbulence 
can introduce an unmodeled dynamic that inconsistently affects the aircraft responses. A frequency- 
varying dynamic uncertainty can be associated with the model to describe nonrepeatible data variations. 

Another source of error between predicted and measured responses is an incorrect assumption of 
flight condition. Flight data sets are sometimes generated at test points that attempt to maintain a constant 
flight condition to match the data sets predicted from a model describing the aeroelastic dynamics at that 
same flight condition. Slight variations in flight conditions while the experimental response is measured 
may cause some discrepancy between the predicted response and the flight data. Parametric uncertainty 
associated with the unsteady aerodynamic model can be used to account for these errors because flight 
condition variations only affect the aerodynamic model and not the structural model. 

The model may accurately represent the mode shapes of the aircraft but have a poor representation of 
the sensor locations. The responses measured by sensors are inherently dependent on sensor location, 
with respect to mode shapes, to determine the magnitude and phase of the signal. Additional errors in 
magnitude and phase are introduced when considering transfer-function estimates generated by signals 
that violate assumptions used in computational algorithms. A frequency-varying dynamic uncertainty can 
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be associated with the output of the plant model to describe errors in both mode shape and sensor 
location predictions. 

The choice of signal-processing algorithms can also introduce errors between predicted and measured 
responses. The Fourier transform, which is the traditional tool for signal processing, assumes severa 
characteristics of the data that may be violated with transient-response aeroelastic data. Filtering and 
wavelet-based algorithms may be used to reduce errors introduced by signal processing and reduce con- 
servatism in the resulting stability margins . 45 * 46 Parameter uncertainty associated with the modal param- 
eters of the linear model may be used to describe some errors in the natural frequencies and dampings 
observed using flight data analyzed by incorrect algorithms. Dynamic uncertainty may also be required to 
describe errors introduced by leakage and aliasing effects. 
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CHAPTER 6 

MODEL VALIDATION OF UNCERTAINTY 


6.1 Model Validation Using the Structured Singular Value 

Robust stability analysis considers the stability of a system subject to a set of modeling errors and 
perturbations represented by a norm-bounded A. A logical question that arises is how to reasonably deter- 
mine this set. This issue is important for computing robust flutter margins because (I can be overly con- 
servative if the uncertainty is excessive beyond the true modeling errors and can be overly optimistic if 
the uncertainty does not sufficiently account for the true modeling errors. 

Model validation algorithms can be used to indicate if an uncertainty description is reasonable with 
respect to flight data. 47 These algorithms consider whether or not a set of data measurements could have 
been generated by a proposed model that includes the nominal dynamics and associated uncertainty oper- 
ators and noise and disturbance signals. 

Uncertainty operators associated with the nominal plant model specifically in the LFT framework can 
be considered by validation algorithms. 48 Frequency-domain algorithms are generated that consider the 
model validation problem in the context of control design. 49 Time-domain approaches are also developed 
that can be solved with convex optimization algorithms for certain uncertainty structures. 50 

A model validation procedure has been formulated that uses a ft, condition to determine if an uncer- 
tainty model is invalidated. 51 This procedure uses frequency-domain transfer- function data to determine 
if some perturbation A € A to the nominal plant could produce the measurements. Consider the block 
diagram for robust stability analysis of systems with measurement y and forcing u signals shown in 
figure 6.1. The model validation question, as applied to uncertainty models, is given in question 6.1.1. 



Figure 6.1. Linear fractional transformation system for robust stability analysis and model validation with 
forcing and measurement signals. 


Question 6.1.1 (uncertainty validation): Is there some frequency-varying A e A for figure 6.1 
such that F U (P, A) could generate the set of observed data y and u? 
n 0 + n z x iij + n w 

Define P(s) e C as a stable transfer-function system matrix such that P e . 

n x n n x n 00 

Partition this matrix into four elements such that P n (s) e C z w , P 22 (s) e C 0 ', and P, 2 (s) and 

p 2 i ( s ) are of appropriate size. P 22 (s) is the nominal-plant transfer function in the absence of uncertainty. 
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Robust stability analysis of this system is determined using the transfer function Pj v which contains 
the feedback relationship between the plant dynamics and the uncertainty operator. The robust stability 
condition given by theorem 3.3.3 requires fi(P j j) < 1. 

The model validation condition uses these elements of P and the finite-energy frequency-domain 
signals y, u e '% 2 where the measurements are scalar functions y, u e C. Formulate the following 

two matrices. 


P 22 = P 22 u -y 

The following theorem formulates the model validation test using |i 51 

Theorem 6.1.2: Given measurements y generated by inputs u, then the system P with associated 
uncertainty set A is not invalidated if 


f*( p ii — P 12 P 22 P 21 ) > * 


(92) 


This condition may seem counterintuitive in that the desired condition for validation is p greater than 
one, while the robust stability condition seeks a value less than one. This condition can be explained by 
considering the following relationship of y = F U (P, A)u as shown in figure 6.1. 

o = [P 22 u-y]+P 21 A(I-P u A) _1 [P 12 u] 

= P 22 + P 21 A(I-P |1 A)- , P I2 

Define the plant P = (P , v p i2> p 2 i’ P 22>- A valu f of |X(P) < 1 implies this system is robustly stable 
to all perturbations A € A. This robust stability of ft(P ) < 1 also implies that the loop gain F u (P, A) is 
not singular for any value A € A. Thus, fi(P ) < 1 would contradict the relationship shown in the above 
equations that requires F u (P, A) to be singular to satisfy the input-to-output relationship of the data. In 
this respect, the model validation test is actually an inverted robust stability test. 

The phrase “model validation” may be misleading. No analytical model can ever be truly validated by 
considering a finite set of experimental data. A model may not be invalidated by the data, but no guaran- 
tee exists that a different data set could not be generated from the physical system that invalidates the 
model. The finite sets of measurement data cannot record the response of the system to an infinite num- 
ber of input signals subject to an infinite number of initial conditions. Theorem 6. 1 .2 reflects this fact by 
explicitly stating the condition only determines if the model is invalidated by the data. The system is 
assumed to be nearly linear, and the data is assumed to sufficiently represent the behavior of the dynam- 
ics, so theorem 6.1.2 represents a model validation test. 



6.2 Validating Norm Bounds for Uncertainty 

The value of p computed using theorem 6.1.2 can be interpreted as a measure of how reasonable the 
uncertainty description is. A value of p = 2 implies the uncertainty can be scaled by 2 before the model is 
invalidated. Alternatively, a value of p = 0.5 implies twice as much uncertainty is required for the model 
to not be invalidated. This interpretation of p, as relating to the size of allowable perturbations, empha- 
sizes the relationship of the model validation condition to a robust stability condition. 

This model validation is used in practice to generate reasonable norm bounds for an uncertainty 
description. The following algorithm can be used to determine a sufficient level of uncertainty required 
such that the model is not invalidated by multiple data sets. A small scalar a > 1 is chosen to scale the 
uncertainty set and increase the amount of allowable errors if the size of A is not sufficient. 

Algorithm 6.2.1 (Model Validation): 

Given frequency-domain data sets {y,, y 2 , .... y n } and (Uj, u 2 , ..., u n }: 

Given frequency-domain transfer-function P with elements P n , P J2 , P 2J , P 22 ; 

Given uncertainty set A with initial norm bound: 

Given update scalar a > 1 : 

valid = FALSE 
while (valid == FALSE) { 
valid = TRUE 
for i = 1 : n { 


P 22 = P 22 u i “ Vi 

if ]l (P 11 -P 12 PJ 2 P 21 XI ( 

valid = FALSE 

} 

} 

P j j = aP j j ( equivalent to A = aA) 

} 

A is the required norm-bounded uncertainty 

The norm bound on the uncertainty set A is not actually increased with this method as denoted by the 
parentheses around the statement A = aA. Scaling the uncertainty and retaining these scalings throughout 
the procedure would be difficult. The algorithm is simplified by always considering the uncertainty is 
scaled such that p < 1 is always the desired result. This scaling of the uncertainty is actually accom- 
plished by the statement P n = aP„, which scales the feedback signals between the plant and uncer- 
tainty operators. The norm bound for the uncertainty that is needed to ensure the uncertainty levels are 
not invalidated by the data is scaled into the plant by the parameter a. Robust flutter margins can be com- 
puted directly from the new scaled plant using algorithm 4.4.6 with a desired p < 1 condition. 
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ai thm f. ? 1 is a straightforward implementation of the model validation test. The p values are 
Algorithm 6.2 1 is a stnugh V invalidated . Mor e sophisticated algorithms could use 

only compared with 1 to determm uncertainty. Also, the values of 1 1 across frequency 

reflect this information. 

a"lpu,e a less conservative uncertainty description that does no. tnvahdate the model. 

Theorem 6 1 2 is only valid for scalar data signals generated by systems with a single input and single 

the number of input and measurement signals, 
damping estimates. 

Another interpretation of algorithm 6.2.1 uses the uncertainty A to bound the magnitude and phase of 
the £Swe Cctions generated by the family of plans F„(P, A). The structured angular value 

« to ex^rimernal data .Lsfer function lies within these analytical bounds a. every frequency. 
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CHAPTER 7 

PROCEDURE FOR THE |I METHOD 


7.1 Model Updating 

Generating a model by analyzing flight data is essential for computing a confident stability analysis. 
A nominal model generated purely from analytical equations of the predicted aircraft dynamics may not 
accurately describe the true aircraft. A model must be generated that accounts for the flight data to ensure 
the predicted dynamics represent the true dynamics. 

The most direct method of generating a model from the flight data is to identify a system model 
entirely from the data measurements. Many system identification algorithms exist that have become 
standard tools for systems and control engineers. 52 Direct application of these methods to aeroelastic sys- 
tems rarely produces an accurate model that accounts for the dynamics of the aircraft. 2 Aeroelastic 
response data is typically of poor quality relative to ground vibration test data because of the low signal- 
to-noise ratio and unobserved dynamics in the response measurements that may drastically lower the 
effectiveness of system identification algorithms. 2 

An alternative method is to use the nominal aircraft dynamical model as an initial estimate to model 
the true aircraft. The flight data are then used to update the elements of this model. Several methods have 
been devised to update an analytical structural model using experimental data. 53 Model updating can be 
performed on the full stress model or a subset computed with Guyan reduction. 54 Generally, considering 
the full model is preferable because the reduction may distribute local errors throughout the entire model 
if an orthogonality condition is violated. 55 

Two basic methods are proposed to update the full structural model using comparisons between 
experimental and predicted data. One method updates the mass and stiffness matrices of the finite- 
element model. 56 This method suffers from lack of physical interpretation of the matrix updates and 
possible numerical conditioning. Another method updates specific parameters in the model. This method 
is accurate for small systems but may require an excessive computational cost for large systems. 

Aeroelastic models have the additional freedom of updating the aerodynamic and the structural 
elements. A method of updating the linear model in a modem control framework has been developed. 32 
This method may be overly conservative for describing nonlinearities, and the corresponding stability 
margins are only accurate for flight conditions near the instability. A parametric identification algorithm 
has been developed that uses flight data to update specific terms in the aerodynamic model through a 
nonlinear optimization. 57 This method suffers with flight data because of unobserved dynamics and the 
low signal-to-noise ratio in the measurements. 

The approach taken in this paper is to update only the uncertainty operators of the robust aeroelastic 
model, using the flight data, and leave the nominal dynamics model unchanged. The model validation 
condition of theorem 6.1.2 is used with the nominal plant P to generate a reasonable uncertainty descrip- 
tion A to associate with P. Figure 7. 1 shows the flow of information through the p method. 

Algorithm 6.2.1 presents a procedure to implement the p method in a manner corresponding to 
figure 7.1. The model update procedure based on algorithm 6.2.1 actually scales the plant but has the 
equivalent effect of scaling the norm bound of the uncertainty description. The algorithm loops over the 
model validation procedure until an uncertainty description is determined that is not invalidated by the 
flight data. The final step is to compute a robust flutter margin from the scaled plant by using a p < 1 con- 
dition as in algorithm 4.4.6. 
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Robust margin 970992 


Figure 7.1. Information flowchart to generate plant and uncertainty operators from a system model and 
flight data with the p method. 


Algorithm 7.1.1 (robust flutter margins with model updating): 

Given nominal plant P: 

Given uncertainty set A associated with P: 

Given input excitation data u: 

Given output response measurement data y: 

Define W > I as the scaling update for A 

while (F (P, A) is invalidated by u, y using algorithm 6.2.1 ) { 

A = WA 

) 

Aflutter is f lutter margin computed from algorithm 4.4.6 

Several advantages exist to using this method as compared to traditional model updating methods. 
The typically poor quality of flight data, in association with aircraft dynamics consisting of many modes, 
makes updating a nominal model difficult. Traditional methods of norm-based update algorithms often 
generate a nonunique set of model updates that have no way to determine which update has the most 
logical physical interpretation. The method of updating the uncertainty operators based on a worst-case 
magnitude avoids this problem. 

Also, this method can work with flight data of varying quality. The updated uncertainty description is 
highly accurate if the data show a high signal-to-noise ratio and much of the dynamics are observed by 
the sensors. If the data do not have these desired characteristics, however, the method can still compute a 
flutter margin. An uncertainty description may be difficult to compute if the data do not indicate the 
aircraft dynamics well, so the model validation procedure will not require a large magnitude for the 
uncertainty operators. The robust model in this situation will closely resemble the nominal model. In this 
way, the p method will always generate a more accurate flutter margin, and at worst, the robust p flutter 
margin will be equivalent to the nominal flutter margin. 
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7.2 Approaches to Use Flight Data 


A flight test generally consists of maneuvers at several test points that may be at identical or different 
flight conditions. The entire flight test program will use many flight tests to measure response data at test 
pom s throughout the flight envelope. The model-updating method that generates uncertainty operators 
can use the entire set of flight data from the different test points. 

Several approaches are formulated to use multiple flight data sets to update the uncertainty descrip- 
lon associated with a nominal plant model. The uncertainty description may be different for each 
approach and the resulting flutter margin will be different for each approach because of the dependence 
of p on the uncertainty set. Two approaches discussed here are denoted as local and global. 

A ! .° Cal a PP roach uses fli § ht data from test points at identical flight conditions. These data are used to 
U !l C , ertaint ^ description for the nominal model at the particular flight condition associated 
with the data. The magnitude of the uncertainty operators is chosen such that a robust model at the single 
flight condition is not invalidated by any of the flight data sets measured at that same flight condition 
with no consideration of data from other flight conditions. 

, :* 0Cal ap P roach sh ,° WS the beneflt ^ independently computing uncertainty descriptions for mod- 
els at different flight conditions. This approach allows each uncertainty description to be more accurate 
because, for example, the flight data may indicate much smaller uncertainty operators are required for 
subsonic plant models even though large uncertainty operators are 

tht model g WOrSt ' CaSC flutter margins wil1 be less conservative because less uncertainty is required for 


Algorithm 7.2.1 outlines the local approach to use flight data and compute flutter margins. 

Algorithm 7.2.1 (local approach): 

Define scalar n f - as the number of flight conditions. 

Define scalar n, as the number of data sets at flight condition i. 

Define vectors |u|, Uj, u^lo* input excitations at flight condition i. 

Define vectors jyj.yf, ■■•,y i ‘| as response measurements at flight condition i. 

Define matrices { P , , P 2 , . . . , P^ } „ plant models at flight condition i. 

for i = 1 : n f { 

Choose initial A ; 
for j= 1 : n - { 

validate F u (P i , A ( ) using yj and uj 
increase size of Aj if necessary to validate 

compute flutter margin 5 f [utter from p(F u (P j , A ( )) 
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A global approach uses the entire set of of the uncer- 

models with ,he ““ ta “ d uncertainty descripti0 " " e 

not invalidated by any of the flight data sets. 

Several advantages and disadvantages exist to ofXeTuheTma^gintecause the 

approach. One disadvantage is a possible arge con dition A single particularly inaccurate plant 

conditions that are better representations of the true dynamics. 

One advantage to this approach, however, ’ J^^J^^t-case flight condition are used 
respect to the entire flight envelope. The conditions. Also, this approach is not very sensitive to 

to generate the uncertainty description o resnonse may only appear in certain data sets. The 

piny measured High, data. A poorly mode <ha. did no, clearly 
local approach would not include n»d- 

observe thrs “ ge e nerated * ilh the glob S a i approach may be more conservahve than the 

S approach, but this approach introduces a corresponding htgher margin of safety. 

Algorithm 7.2.2 outlines the global approach to use flight data and compute flutter margins. 

Algorithm 7.2.2 (global approach): 

Define scalar n f as the number of flight conditions. 

Define scalar n, as the number of data sets at flight condition i. 


Define vectors I u| , U u"‘ 1 as input excitations at flight condition i. 


12 n i| 

Define vectors <y { ,y r r as 


response measurements at flight condition 1 . 
Define matrices {P,.P 2 , .... V as plant models at flight condition i. 


Choose initial A 
for i = 1 : Of { 

for j = 1 •- nj { 

validate F (J (P i , A) using y{ and u- 
increase size of A if necessary to validate 

} 

} 

for i = 1 : n f { 

compute flutter margin §' flutler f rom ^ F u( P i’ 

} 
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Hybrid approaches have also been formulated that mix the local and global approaches. One straight- 
forward hybrid approach is to generate an uncertainty description using all data from a small range of 
flight conditions. This approach can be useful for separately considering sets of plant models that are 
generated using different techniques. For example, the model-generating package used for this paper 
computes all subsonic plant models with a doublet-lattice algorithm and the supersonic models are gener- 
ated with constant-panel algorithms. A hybrid approach could be used to reflect this knowledge and 
consider groups of subsonic, supersonic, and transonic plants independently. 

The approaches outlined here are certainly not exhaustive. A weighted-norm approach can be formu- 
lated that uses flight data from the entire flight envelope but depends heavily on a particular subset of that 
data. Other approaches could concentrate on particular dynamics through modal filtering techniques to 
generate separate uncertainty descriptions for individual modes . 58 
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CHAPTER 8 

ROBUST FLUTTER MARGINS OF THE F/A-18 SYSTEMS RESEARCH AIRCRA 
8.1 Aircraft Flight Test 

Mojave desert near Edwards, California. 



Figure 8.1 F/A-18 systems research aircraft. 


Flutter testing was initiated on the SRA because of a major left-wing structural modification that 
allows testing o/several hydraulic and electromechanical aileron actuator concepts. The increased size 
and weight of these actuators required large and heavy items to be installed, replacing a fitting called a 
^LTh ffthlt supports the aileron hinge, actuator, and a fairing. These stmctural modifications changed 
stiffness and dampfng properties and added approximately 20 lbm to the wing. Dependency of the aileron 
aeroelastic behavior on actuator dynamics warranted the flutter tests. 

The results presented in this paper are associated with the F/A-18 SRA with the left-wing structural 
modification More than 30 flight tests were conducted in two sessions, one session between ep em 
m* and Febmt^ 1995 and one session between June 1995 and July 1995. A total of 260 data sets are 
measured from test points at various conditions throughout the flight envelope. 

The flight flutter test program used a wingtip excitation system developed by Dynamic Engineering 
Incorporated (DEI). This exciter is a modification of an excitation system used for F- 
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research. 60 The system consists of a wingtip exciter, an avionics box mounted in the instrumentation bay, 
and a cockpit controller. 

Aerodynamic forces are generated by the wingtip exciter. This exciter consists of a small, fixed, aero- 
dynamic vane forward of a rotating, slotted, hollow cylinder. Rotating the cylinder varies the pressure 
distribution on the vane and results in a wingtip force changing at two times the cylinder rotation fre- 
quency. The magnitude of the resulting force is determined by the amount of opening in the slot. 
Figure 8.2 shows the F/A-18 aircraft with a left-side wingtip exciter in the aft position. 



Figure 8.2 Left wing of the F/A-18 systems research aircraft with wingtip excitation system. 


The cockpit controller commands a frequency range, duration, and magnitude for the wingtip excita- 
tion signal. Frequency-varying excitation is generated by changing the cylinder rotation frequency with 
sine sweeps. Each wingtip exciter is commanded to act in phase (0°) or out of phase (180°) with each 
other. Ideally, the in-phase data excite the symmetric modes of the aircraft, and the out-of-phase data 
excite the antisymmetric modes. In reality, nonideal phasing of the exciters often caused a mixture of 
symmetric and antisymmetric modes to be excited. 

Flight data sets were recorded by activating the exciter system at a given flight condition. The aircraft 
attempted to remain at the flight condition throughout the series of sine sweeps desired by the controller. 
The sine sweeps were restricted to within 3 Hz and 35 Hz, and smaller ranges were sometimes used to 
concentrate on a specific set of modal responses. Multiple sets of either linear or logarithmic sweeps were 
used with the sweep frequency increasing or decreasing. Table 8.1 shows the flight test matrix of excita- 
tion parameters and flight conditions for the F/A-18 SRA flutter program. 

Accelerometers are available at several points on the aircraft to record modal responses. Each wing 
has a sensor on the aileron and at the forward and aft position on the wingtip. Additional accelerometers 
are located on each vertical tail and horizontal stabilator. 
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Table 8.1. Operating and flight conditions for the F/A-18 SRA flight flutter test. 


Mach 

0.54,0.65,0.70,0.80,0.85,0.90,0.95, 1.05, 1.2, 1.4, 1.6 

Altitude, ft 

10K, 30K, and 40K 

Exciter configuration 

Symmetric, antisymmetric, and independent 

Exciter position 

Both forward, both aft, and left aft/right forward 

Exciter force level 

Low and high 

Sweep type 

Linear and logarithmic 

Sweep duration, sec 

15, 30, and 60 

Sweep range, Hz 

3-35, 3-12, 3-25, 25-35, 35-3, 3^0 

Sweep multiples 

1, 2, and 4 sweeps for each maneuver 


8.2 Aeroelastic Data Analysis 

Aeroelastic flight data generated with the exciter system are analyzed by generating transfer functions 
from the excitation force to the sensor measurements. These transfer functions are generated using stan- 
dard Fourier transform algorithms. Several inherent assumptions associated with Fourier analysis exist 
that are violated with the flight data. The assumptions of time-invariant stationary data composed of sums 
of infinite sinusoids are not met by this transient-response data. ’ The amount of uncertainty required 
to validate the model can be reduced by using wavelets to accurately process the data but are not consid- 
ered in this paper. 

Figure 8.3 shows an example of dynamical variations that are observed during flight testing. Two 
transfer functions generated from different data sets clearly show some variation. The natural frequency 



Figure 8.3. Flight data transfer functions from symmetric excitation to left wing forward acceler- 
ometer for Mach = 0.8 and 30,000 ft demonstrating variation in modal frequency and damping. 
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and modal response are different as observed by these two data sets, although each test point used the 
same sweep procedure at identical flight conditions of Mach 0.8 and an altitude of 30,000 ft. One 
possibility to account for such variations is the change in mass, a maximum of 10 percent, between the 
heavyweight and lightweight conditions caused by fuel consumption. Another possibility is the nonideal 
phasing between the exciters causing excitation that is not purely symmetric. 

The deviations between the modal responses in figure 8.3 are used to generate uncertainty for the 
symmetric wing first-bending mode. Similar deviations observed in modal responses are used to generate 
uncertainty for each mode. Uncertainty levels in modes that are not observed are assumed to be similar to 
those of the observed modes. 

Figure 8.3 shows modeling uncertainty caused by variations between the analytical model and the 
aircraft dynamics that change during a flight test. Additional uncertainty is introduced because of 
improper operation of the exciters. Several behaviors are observed because of these improper operations. 

The most noticeable of these improper behaviors occurs as a function of dynamic pressure. The 
motors and rotating cylinder assemblies bind at high dynamic pressures, causing erratic phase relation- 
ships between the individual exciters on opposite wingtips. Figure 8.4 shows a representative plot of the 
phase difference in degrees between the left and right exciters operating in symmetric mode for a sweep 
taken at a test point with q = 356 lbf/ft 2 . This phase plot, which does not maintain the desired 0° phase 
difference but is relatively constant and well-behaved, is in direct contrast to the phase difference at a test 
point of q = 825 lbf/ft 2 (fig. 8.5). 



Figure 8.4. Phase difference in degrees between left and right exciters for sweep range 3-30 Hz at 
Mach = 0.90 and 30,000 ft at q = 356 lb/ft 2 . S 
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Figure 8.5. Phase difference in degrees between left and right exciters for sweep range 3 30 Hz at 
Mach = 0.90 and 10,000 ft at q = 825 lb/ft . 


Another improper behav.or is the nonrepeatibility observed w.th the data because of tnmw m 
transfer functions using different data sets recorded at identical flight condmons. figures 1 8 ’ 
show an example of such a variation observed for frequency-increasing and -decreasing sweeps. The p 
method can use the data despite these nonrepeatabilities that would adversely affect traditional system 

identification approaches. 

The data analysis process is further degraded by a poor measurement of the excitation force used to 
generate the modal responses. The excitation force is not directly measured, but a strain-gage measure- 
ment is used to approximate this force. The strain gage records lateral shear strain at the exciter vane root 
This measurement is assumed to be representative of the vertical shear and spanwise-moment load at the 
wingtip rib. Time-frequency analysis using wavelets clearly shows this strain-gage signal to have some 
harmonic resonances from the structural response. 

The effect of the poor approximation of input force and the erratic behavior of the exciters is to 
reduce the quality of the flight data. Methods relying on system identification fail to accurately use t e 
data to predict a flutter boundary. 57 The p method is able to account for the data anomalies by including 

greater levels of uncertainty. 
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Figure 8.6. Flight data transfer function from symmetric excitation to left wing forward accelerometer for 
Mach = 0.8 and 30,000 ft with frequency increasing sweep from 3-35 Hz. 



Figure 8.7. Flight data transfer function from symmetric excitation to left wing forward accelerometer for 
Mach = 0.8 and 30,000 ft with frequency decreasing sweep from 35-3 Hz. 
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8.3 Analytical Model 

A finite-element model of the SRA is used to compute the aeroelastic dynamics and modal 
characteristics of the aircraft. The generalized equations of motion are used to derive a linear, finite- 
dimensional state-space model of the aircraft. This model contains 14 symmetric structural modes, 
14 antisymmetric structural modes, and 6 rigid-body dynamic modes. The control surfaces are not active, 
and no control modes are included in the model. Table 8.2 shows natural frequencies of the elastic struc- 
tural modes. These modal frequencies are computed for the aircraft with no aerodynamics considered. 


Table 8.2. Elastic natural frequencies of the structural finite-element model of the fully 
fueled F/A-18 SRA with left-wing modifications and no wingtip excitation system. 


Mode 

Symmetric, Hz 

Antisymmetric, Hz 

Wing first bending 

5.59 

8.84 

Fuselage first bending 

9.30 

8.15 

Stabilator first bending 

13.21 

12.98 

Wing first torsion 

13.98 

14.85 

Vertical tail first bending 

16.83 

15.61 

Wing second bending 

16.95 

16.79 

Wing outboard torsion 

17.22 

— 

Fuselage second bending 

19.81 

18.62 

Trailing-edge flap rotation 

23.70 

23.47 

Fuselage torsion 

- 

24.19 

Launcher rail lateral 

- 

24.35 

Stabilator fore and aft 

28.31 

28.58 

Wing second torsion 

29.88 

29.93 

Aileron rotation 

33.44 

- 

Aft fuselage torsion 

- 

37.80 

Aileron torsion 

38.60 

— 

Wing pitch 

- 

39.18 

Wing third bending 

43.17 

— 


The predicted flutter results for this aircraft are computed from the finite-element model using the p-k 
method. A detailed explanation of the F/A-18 SRA modeling process and flutter analysis using traditional 
methods has previously been published. 

The doublet-lattice and constant-panel methods are used to compute the frequency-varying unsteady 
aerodynamic forces for several subsonic, transonic, and supersonic Mach numbers. Force matrices for 
Mach 0.80, 0.90, 0.95, 1.10, 1.20, 1.40, and 1.60, are available. The unsteady aerodynamic forces are 
computed as a function of reduced frequency, k. 
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(94) 


k = 



The reduced frequency is a function of the true frequency (to), the true velocity (V), and the mean aero- 
dynamic chord (c). Aerodynamic forces generated for ten reduced-frequency points between k = 0.0001 
and k = 4.0000 are sufficient for flutter margin computation for this aircraft. Figure 8.8 shows the magni- 
tude and phase of a sample calculation of the unsteady aerodynamic forces as a function of reduced 
frequency. An additional ten points computed between k = 0.01 and k = 1.00 demonstrate the function is 
smooth across frequency, so the original set of ten points is sufficient to describe the behavior of unsteady 
aerodynamic forces. 


Log 10 o 
magnitude 


Ten points between 
reduced frequency 
K = 0.0001 and K = 4 
K = 0.01 and K = 1 




Figure 8.8. Magnitude and phase of the frequency-varying unsteady aerodynamic forces. 


The unsteady aerodynamic forces are fit to a finite-dimensional state-space system. The system 
identification algorithm is a frequency-domain curve-fitting algorithm based on a least-squares minimi- 
zation. A separate system is identified for each column of the unsteady forces transfer-function matrix. 
Fourth-order state-space systems are used for each column of the symmetric forces, and second-order 
state-space systems are used for each column of the antisymmetric forces. These systems are combined to 
form a single multiple-input and multiple-output state-space model of the unsteady aerodynamics forces, 
previously designated Q(s). This model has 56 states for the symmetric modes and 28 states for the anti- 
symmetric modes. 
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The analytical aeroelastic model has inputs for symmetric and antisymmetric excitation forces. The 
excitation force is assumed to be purely symmetric or antisymmetric. Ten outputs exist that correspond to 

the sensor locations on the aircraft. 


8.4 Uncertainty Description 


Noise and uncertainty operators are introduced to the linear aeroelastic model to account for varia- 
tions between the analytical model and the actual aircraft. Standard analysis of the linear model is used to 
determine the framework for how uncertainty operators enter the system. Two uncertainty operators an 
a single noise input are used to describe the modeling uncertainty m the linear aeroelastic model The 
magnitude of each uncertainty operator and the noise level is determined from both reasoning of the mod- 
eling process and analysis of the flight data. 

An uncertainty operator, A a , is associated with the state matrix of the F/A-18 linear model This 
uncertainty models variations in the natural frequency and damping values for each mode. State m 
uncertainty can account for errors in coefficients of the equations of motion and changes in the aircraft 
dynamics caused by parameter variations such as mass consumption during flight. 


The A» is a structured, diagonal matrix with real, scalar parameters as elements. Separate elements 
are used to affect each modal response and time lag in the state matrix. The modal response uncertainty 
parameters are each repeated two times, and each time-lag uncertainty appears once on the diagonal. 

Each repeated modal uncertainty parameter affects natural frequency and damping by allowing varia- 
tion in the state matrix elements. Consider formulating the state matrix as a block diagonal with a 2-by-2 
block representing each mode. The real component of a modal eigenvalue, e r , is the diagonal component 
of each block, and the imaginary part, ep is arranged on the off-diagonal positions. Define \ as the block 
determining the natural frequency, w^, and damping, Cj, of the i mode. 



- 

/ 2 2 

_ 

e r e . 

^ w i = y e r +e i 


_- e i e r_ 

Ci = -e r /», 


(95) 


Scalar weightings, w r and Wj , are used to affect the amount of uncertainty in each matrix element. 
The amount of variation in the matrix elements and, correspondingly, the amount of variation in thenatu- 
ral frequency and damping, are determined by the magnitude of these scalar weightings. Define e r and 
^ as elements of the state matrix affected by an uncertainty parameter 6. 


e-, = e r ( 1 ± w, 5) 

Aeroelastic modes typically show low damping values caused by the real component being smaller 
than the imaginary component. Because linear modeling techniques often identify the natural frequency 
better than the damping value, the weighting for the real component is expected to be larger than that for 
the imaginary component. This weighting is reflected by the observed modal parameters in the flight data. 
The natural frequencies show variations of ±5 percent from the theoretical model, but the uncertainty in 
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the damping needs variations of approximately ±15 percent to validate all the flight data. The scalar 
weightings are chosen accordingly. 


w r = 0.15 
Wj = 0.05 


(97) 


The flight data are only able to determine uncertainty levels for the modal parameters of the experi- 
mentally observed modes. The uncertainty levels in the unobserved modes are assumed to be consistent 
with these values. Parametric uncertainty is introduced for each modal block in the state matrix, affecting 
observed and unobserved modes, using the weighting values given above. 


The single scalar blocks of A a are normalized to one by weighting their effect on the time lags of the 
state matrix. Variations observed in the flight data are used to determine that a weighting of w, = 0. 15 
is required to admit 15-percent variation in the time lags. ag 

The second uncertainty operator, A in , is a multiplicative uncertainty on the force input to the linear 
model. This uncertainty is used to cover nonlinearities and unmodeled dynamics. The linear model 
contains no dynamics at greater than 40 Hz, so the high-frequency component of this operator will reflect 
this uncertainty. This operator is also used to model the excitation uncertainty caused by the exciter sys- 
tem. Analysis of the flight data indicates the input excitation signals rarely had the desired magnitude and 
phase characteristics that they were designed to achieve. The low-frequency component of the input 
uncertainty reflects the uncertainty associated with the excitation system used to generate the flight data. 
The frequency- varying transfer function for weighting the input uncertainty is given as Wj . 


W = 5 s + 

,n s + 5000 


(98) 


A noise signal is included with the sensor measurements. Knowledge of the aircraft sensors is used to 
determine that a level of 10-percent noise is possible in the measured flight data. An additional noise may 
be included on the force input caused by the excitation system, but the input multiplicative uncertainty is 

sufficient to describe this noise. Figure 8.9 shows the block diagram for the aeroelastic model with the 
uncertainty operators. 



Figure 8.9. F/A-18 robust stability block diagram. 
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The flight data used to validate this uncertainty structure cover a large range of flight points. The 
entire set of 260 flight maneuvers throughout the flight envelope is considered. Using a single ^certain y 
descripflon over the entire flight envelope may be conservative. Certainly assuming the linear models are 
more accurate at subsonic and supersonic speeds than at transonic speeds because of flow irregularities 
earned by shock formation is reasonable. Additionally, the flight data from the exciter system should be 
at lubsorhc speeds than at supersonic speeds because of the degraded performance at high dynamic 
the analysis process is simplified by considering a single set 
tors This process is equivalent to formulating the worst-case uncertainty levels at the wors -case flight 
condition and assuming that amount of uncertainty is possible for the remaining flight condition . 

8.5 Nominal and Robust Flutter Pressures 

Flutter pressures are computed for a linear model with the associated modeling uncertainty structure 
using the tt analysis method.® Linear systems for symmetric and antisymmetric structure modes are sep- 
arated for ease of analysis. These systems can easily be combined and analyzed as a single sys em, ow 
ever eigenvector analysis would be required to distinguish which critical flutter modes are symme n 
and whfch are antisynLetric. Each system contains the same number of structural modes (14), and 
uncertainty descriptions are identical for each linear model. 

The system shown in figure 8.9 contains three uncertainty blocks. The parametric uncertainty cover- 
ing variations caused by dynamic pressure, 5^, is a scalar parameter repeated 14 times once or each 
elastic mode. The parametric uncertainty block affecting the modal parameters, A a , is a diagonal matrix 
with dimension equal to the number of states. Separate scalars along the diagonal represent uncertainty in 
each elastic mode, each mode in the aerodynamic force approximation, and each lag term. The uncer- 
tainty parameters for the modes are repeated two times; the parameters for the lag terms are single 
scalars Define 5- as the i th uncertainty parameter for the system using n m modes and n, lag terms. The 
input-multiplicative uncertainty block, A in , is a scalar for this single-input plant model because symmet- 
ric excitation is analyzed separately from antisymmetric excitation. 

The parametric uncertainty parameters represent changes in elements of the state-space model. The 
variation of 8- between ±1 admits dynamic pressures between 0 < q < 2q 0 - Allowing the modal uncer- 
tainty parameters, 5 p 5 n , to vary between ±1 allows 5-percent variation in the imaginary part of the 
natural frequency and 15-percent variation in the real part. This variation corresponds to approximately 
5-percent variation in the natural frequency and 15-percent variation in the damping value of each mode. 
These parameters are real quantities. The input-multiplicative uncertainty contains magnitude and phase 
information and is treated as a complex linear time-invariant uncertainty. 

Nominal flutter pressures are initially computed by ignoring the modal and input uncertainties. The 
structured singular value is computed only with respect to the parametric uncertainty, allowing a range of 
dynamic pressures to be considered. Robust flutter pressures are computed with respect to the structure 
uncertainty set. A, described above using the structured singular value. Table 8.3 shows tradmonal flutter 
pressures computed using the p-k method and nominal and robust flutter pressures computed with \x. 
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Table 8.3 Nominal and robust flutter pressures computed with p-k and p methods. 


Mach 

Symmetric, lbf/ft 2 

Antisymmetric, lbf/ft 2 

^p-k 

^ nom 

^ rob 

Vk 

^ nom 

^ rob 



3168 



4593 

3648 



2706 

2575 


3057 

2944 


2430 

2388 

2329 


2751 

2572 



5676 



3265 

2827 


2469 

2454 

2327 


2893 

2653 

1.40 

3528 

3432 



4439 

4191 



4487 

3996 


5870 

4536 


The nominal flutter dynamic pressures computed using the |i method can be directly compared with 
those computed using the traditional p-k method. Each of these flutter solutions is based on an analytical 
model that has no consideration of modeling uncertainty. 

The nominal flutter points for the symmetric modes match closely with the p-k method throughout 
the flight envelope. The subsonic and supersonic cases show an especially good correlation with the p-k 
flutter points. For each of these flight regions, the (I analysis flutter dynamic pressures are nearly identical 
(within 1 percent) to the p-k method flutter dynamic pressures. The transonic case at Mach 1.1, however, 
shows a slight difference between the two methods. The p method computes a flutter point that is greater 
than the p-k method. In each Mach regime (subsonic, supersonic, or transonic), the nominal flutter points 
are within 5 percent for the two methods. 

The antisymmetric modes show a similar relationship between the flutter pressures computed with 
the p and p-k methods. The subsonic and supersonic flutter points are within 5 percent for the two 
methods, but a greater deviation exists at the transonic condition. The p computes a flutter pressure at 
Mach 1.1 that is 40 percent lower than the p-k method indicates. 

The nominal flutter points for the p and p-k methods show the greatest difference for both the 
symmetric and antisymmetric modes at the transonic case. The aerodynamics at Mach 1 . 1 are more 
difficult to accurately model than at either subsonic or supersonic speeds. Numerical sensitivity to repre- 
sentations of the unsteady aerodynamic forces causes differences in the nominal flutter pressure. 

The robust flutter pressures computed using the p method have lower dynamic pressures than the 
nominal pressures, which indicates the expected conservative nature of the robust computation. These 
new flutter points are worst-case values for the entire range of allowed uncertainty. The subsonic and 
supersonic flutter boundaries are not greatly affected by the uncertainty set. In each of these cases, the 
robust flutter point is within 10 percent of the nominal flutter point. 

The flutter boundary at the transonic case, Mach 1.1, demonstrates significant sensitivity to the 
modeling uncertainty. The robust flutter dynamic pressures are approximately 70 percent of the nominal 
flutter pressures. This difference is explained by considering the rapid transition of critical flutter bound- 
aries near this region. The critical flutter frequencies and the flutter dynamic pressure widely vary 
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between Mach numbers slightly lower and higher than transonic_The small ' 

tainty is enough to cause the worst-case flutter mechanism to shift, and the plant assumes characte 
more consistent with a nontransonic regime. 

Table 8 4 shows the modal natural frequencies for the critical flutter modes^ The frequencies 
enmnuted usina the p-k method and the p analysis method are close throughout the flight envelope for 
Ld antisymmetric modes Frequencies for the robust flutter solutions are slightly 
different than the nominal flutter frequencies because of the modeling uncertainty that allowed 5-percent 
variation in the modal natural frequencies. 


Table 8.4. Frequencies of unstable modes 
computed with p-k and (i. methods. 

associated with nominal and robust flutter pressures 



Symmetric, Hz 


Antisymmetric, Hz 


Mach 

M 

1 

o. 

3 

®nom 

®rob 

® P -k 

®nom 

®rob 

0.80 

8.2 

7.6 

7.7 

9.0 

9.1 

9.1 

0.90 

7.4 

7.3 

7.3 

9.2 

9.1 

9.2 

0.95 

6.8 

6.9 

6.9 

9.1 

9.2 

9.2 

1.10 

12.1 

13.2 

13.0 

28.6 

28.0 

28.3 

1.20 

26.5 

21 A 

27.4 

26.9 

28.9 

28.9 

1.40 

28.1 

28.1 

28.1 

30.4 

31.7 

31.8 

1.60 

28.9 

30.1 

30.1 

32.8 

32.3 

32.1 


Table 8 5 shows subcritical flutter pressures computed with the p and p-k methods. Only nominal sub- 
critical pressures are detected with fl, because the robust pressures are always worst-case critical pressures. 


Table 8.5. 

Subcritical nominal flutter pressures computed with p-k and p 

methods. 


Symmetric, lbf/ft 2 

Antisymmetric, lbf/ft 

Mach 

^p-k 

^ nom 

*7 p-k 

^ nom 




4700 

4583 

0.90 



5300 

5093 

0.95 

7450 

6919 



1.10 



6050 

6001 

1.20 

5400 

5003 

8400 

7943 

1.40 

8970 

8959 



1.60 

8400 

8843 




The it analysis computes subcritical flutter pressures within 10 percent of the p-k method for both the 
symmetric and antisymmetric modes. The p method is even able to detect the subcritical flutter hump 
mode occurring for antisymmetric excitation at Mach 0.9. 
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8.6 Nominal and Robust Flutter Margins for a Flight Envelope 

The nominal and robust flutter pressures shown in table 8.3 are used to compute flutter margins for a 
flight envelope. These flutter margins are determined by the distance between the flutter pressure and the 
dynamic pressure at the edge of the flight envelope at a particular Mach number. 

Table 8.6 shows the flight conditions for several points at greater than Mach 0.8, which have the high- 
est dynamic pressures in the defined flight envelope and the nominal and robust flutter margins. All of 
these flutter margins have positive values, indicating the envelope is free of flutter instabilities. The criti- 
cal flutter pressures are simply the sum of the flight envelope dynamic pressures and the flutter margins. 


Table 8.6. Nominal and robust flutter margins for a flight envelope expressed as allowable 
increase in dynamic pressure computed with p-k and p methods for symmetric and 
antisymmetric modes. 


Flight envelope 

Symmetric, lbf/ft 2 

Antisymmetric, lbf/ft 2 

Mach 

q. 

r P- k 

r 

nom 

T rob 

r P- k 

r 

nom 

r rob 

0.80 

948 

2412 

2220 

1961 

3652 

3645 

2700 

0.90 

1200 

1500 

1506 

1375 

1950 

1857 

1744 

0.95 

1337 

1093 

1051 

992 

1263 

1414 

1235 

1.10 

1491 

3909 

4185 

2629 

4009 

1774 

1336 

1.20 

1438 

1031 

1016 

889 

1412 

1455 

1215 

1.40 

1252 

2276 

2180 

1782 

3348 

3187 

2939 

1.60 

1027 

3443 

3460 

2969 

4673 

4843 

3509 


The flutter margins are often expressed in terms of the percent change, IT, in knots of equivalent 
airspeed, KEAS, between the flight envelope and the flight condition corresponding to the critical flutter 
pressure. Table 8.7 shows these percentage values. 


Table 8.7. Nominal and robust flutter margins for a flight envelope expressed as percent of allow- 
able increase in KEAS computed with p-k and p methods for symmetric and antisymmetric modes. 


Flight envelope 


Symmetric 


Antisymmetric 

Mach 

KEAS 

V 

^nom 

n rob 

V 

^ nom 

V 

0.80 

529 

88 

82 

75 

120 

120 

96 

0.90 

595 

50 

50 

46 

62 

59 

56 

0.95 

628 

35 

33 

31 

39 

43 

38 

1.10 

663 

90 

95 

66 

92 

48 

37 

1.20 

652 

31 

31 

26 

41 

41 

35 

1.40 

608 

67 

65 

55 

91 

88 

82 

1.60 

550 

108 

109 

97 

135 

139 

109 
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A flutter margin yielding a percentage II = 15 for allowable increase in KEAS is commonly used to 
. , fi; CT u t pnvplone The flutter margins shown in table 8.7 demonstrate the indicated 
n!^nvc“« f/A- 18 SRA is acceptable because the nominal and robust nutter margins are all 

c eare ^Ts percent for symmetric and antisymmetric modes. The smallest flutter ™rg.n occurs for 
C— lie dynamics a. Mach 1.2, which presents a robust flutter margtn of only 26 percent, 
however, even this condition is greater than the desired 15-percent margin. 

The dvnamic pressures at which flutter occurs are converted into altitudes, commonly known as 
matlhed-point solutions, using standard atmospheric equations. Figure 8.10 

for the symmetric modes and figure 8.11 shows them for the antisymmetric modes. The flight P 

“/aTs^ is shownl these plots, along with the required 15-percent flutter boundary for 

military aircraft. 

Figures 8. 10 and 8. 1 1 show several short solid lines to indicate the p-k flutter solutions J hr ^ hou ‘ * e 
flight regime Each of these short solid lines represents the flutter points caused by a specific mode_ 
tefpoinfs for the symmetric modes shown in figure 8.10 show four solid lines indicating three different 
critical flutter modes for the considered range of Mach numbers and a subcntical flutter mode occurring 
at supersonic Mach numbers. In figure 8.11, the antisymmetric modes show Ihe onset ° f ““ 
threedifferent critical modes and three subcritical flutter modes throughout the flight envelope. Table 
shows the frequencies of the critical flutter modes. 

The subsonic flutter altitudes for symmetric and antisymmetric modes < ‘ ^ 95 ^ 
teristic The nominal flutter boundary shows a significant variation from Mach 0 80 to Mach 0.95 that is 
caused by sensMvity to Mach number for the dynamics associated with the critical flutter mode. Th 
robust flutter boundary indicates the sensitivity of the plant to errors and the worst-case perturbation The 
higher altitude for the nominal flutter boundary at Mach 0.81 than for Mach 0.80 is reflecte in ® ^ 
conservatism associated with the robust flutter boundary. Similarly, slight variation of Mach number 
Mach 0.95 is not expected to decrease the nominal flutter margin, so less conservatism is associated with 

the robust flutter boundary. 

An interesting trend is noticeable for the symmetric mode robust flutter points in figure 8.10 at the 
supersonic Mach numbers. The flutter mechanism results from the same modes from Mach 1.2 
Mach 1.6 with some increase in frequency. Similarly, the altitudes of the nomm utter ^^rg^s s ow 
little change for these Mach numbers. The aeroelastic dynamics associated with the critical flutter mode 
are relatively unaffected by the variation of Mach over this range, and consequently, each flutter bound- 
ary has the same sensitivity to modeling errors. 

The robust flutter margins for the antisymmetric modes at supersonic Mach numbers show a slightly 
different behavior than the symmetric mode flutter margins. The flutter instability is again caused by a 
single mechanism from Mach 1 .2 to Mach 1 .6, with similar frequency variation for the symmetric modes. 
The robust flutter margins demonstrate a similar sensitivity to modeling errors at Mach 1 2 land l Ma ich 
but at Mach 1.6, a greater sensitivity is shown. The greater conservatism at Mach 16 may indicate 
impending transition in flutter mechanism to the subcritical mode at slightly higher Mach number. 

Figures 8 10 and 8.1 1 show dark solid lines representing the required boundary for flutter P 01 ^ ts _ A11 
nominal and robust flutter points lie outside this region, indicating the flight envelope should be safe from 
flutter instabilities. The robust flutter boundaries computed with p indicate more danger of encoun ering 
flutter exists than was previously estimated with the p-k method. In particular, the robust flutter margin for 
symmetric excitation at Mach 1.2 lies considerably closer to the boundary than the p-k method indicates. 
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Nominal p-k margin 

Robust (i margin 

o Nominal ji margin 


Mach 


971000 


Figure 8.10. Nominal and robust flutter points for symmetric modes. 



8.7 Computational Analysis 

The p analysis method of computing flutter margins presents significant analytical advantages 
because of the robustness of the resulting flutter margin, but the method also has several computational 
advantages over the p-k method. The p algorithm requires a single, linear, aeroelastic, plant model at a 
given Mach number to compute critical and subcritical flutter margins. An entire set of flutter margins 
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can be easily generated in a few minutes using a standard engineering workstation and commercially 
available software packages. 

The p-k method is an iterative procedure that requires finding a matched-point solution. The aircraft 
is analyzed at a particular Mach number and air density. The airspeed at these conditions resulting in a 
flutter Instability^ is computed. This airspeed, however, often does not correspond to the unique airspeed 
determined by that Mach number and air density for a standard atmosphere. Various air densities are used 
to compute flutter solutions, and the corresponding airspeeds are plotted. Figure 8.12 shows an examp 
of an airspeed plot for flutter. 

Figure 8 12 shows vertical lines representing two antisymmetric modes that may flutter at Mach 1.4. 
The p-k method computes a flutter solution at the airspeed indicated where the modal line crosses the 
stLard-Ttmosphere curve. This flutter solution is difficult to compute from only a few air density com- 
putations Typically, several air densities are used to compute airspeed flutter solutions and a line is 
extrapolated between the points to determine the matched-point solution at the standard-atmosphere 
crosXg point The nonlinear behavior shown for mode 1 near the standard-atmosphere crossing point 
indicates an accurate flutter boundary would be extremely hard to predict unless many solutions are com- 
puted near the true matched-point solution. 

The p-k method can also have difficulty predicting the subcritical flutter margins. The second mode 
in figure 8.12 may or may not intersect the standard-atmosphere curve. More computational analysis is 
required to determine the behavior of this mode at high airspeeds. The \i analysis method accurately 
detects the critical and subcritical flutter margins without requiring expensive iterations. 
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CHAPTER 9 

RESEARCH EXTENSIONS 


9.1 Sensitivity Analysis 

Stability margins, especially for high-performance aircraft, may be extremely sensitive to the aircraft 
dynamics and flight conditions. Analytical models must be accurate to ensure the estimated stability 
properties represent the true aircraft stability properties. Such accuracy is difficult to achieve for complex 
modern aircraft designs because components such as flexible structural elements and flight regimes such 
as lgh angles of attack are complex to model. The sensitivity of the stability margins to changes in the 
mo lei indicates some degree of confidence in the predicted aircraft stability properties. Margins with low 

sensitivity are relatively trustworthy; margins with high sensitivity may be inaccurate in the presence of 
modeling errors. F 

Sensitivity of flutter margins with respect to aeroelastic dynamics is traditionally considered using 
statistical approaches. Elements of an analytical model are subjected to random perturbations with known 
s atistical properties, and changes in stability margins are monitored. A state-space model is used to 
compute sensitivity to modal variations with a first-order approach . 12 A finite-element model is used to 
estimate sensitivity to higher-order perturbations . 16 These methods are based on Monte Carlo simulations 

at can be extremely time-consuming and provide no guarantees of sensitivity to perturbations not 
explicitly used in the simulation. * H n!> not 

The n method can be directly applied to compute a measure of the sensitivity of the stability proper- 
les to variations in the analytical model. The conservatism introduced by the worst-case stability margin 
is exactly a measure of the sensitivity to variations on the order of the uncertainty operators »!y 
een nominal and robust margins indicates low sensitivity with respect to perturbations included in A 

Etefmite advantages exist to computing sensitivity with |n instead of statistical approaches The 
exceedingly high number of perturbations needed for Monte Carlo simulations of complex systems is 

AkoTh ^ \ Smg e calculatlon because 0 P er ator theory ensures ft considers all perturbations A € A 
Also, the sensitivity is guaranteed to be worst-case with respect to these perturbations. 

9.2 On-Line Stability Tracking 

Flight flutter testing incurs dramatic time and costs for safely expanding the flight envelope to ensure 
lnstabll dies are encountered. Traditional flight test methods for determining flutter margins 
track estimates of modal damping obtained from flight data are inefficient and dangerous 2 Muldple 
ata sets must be taken at each test point to ensure sufficient excitation of critical modal dynamics Also 
the envelope must be expanded in small increments because damping only indicates stability at the 

onTne ITth °T ^ r? rdiably extra P olated to consider other flight conditions. Several proposed 

on-line methods use modal filtering, parameter identification, and envelope analysis to increase confi 

These rae,hods ’ " ke dampi " 8 ~ 


The J l method is extended to compute stability estimates on-line during a flight test The procedure 
uses flight data from the current test point to update the uncertainty description and compute a new flutter 
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margin 67 Local and global approaches are formulated to reflect different degrees of confidence rn the 
model and flight data. 

The concept of a flutterometer is introduced as a flight aabd!" 

condition. 38 The on-line application of H is a js ^^d to greatly improve the effi- 
predictor and not merely a sta “'^ , can be red uced because the analysis time at each test 


9.3 Robust Aeroservoelastic Stability Analysis 


Aeroservoelastic dynamics couple structural 

Correct m^gi^of^roM^oelastic^tabUU^may result from analysis .ha, does no, account for these 
. . 69 

model uncertainties. 

Several traditional approaches are used to a^alyzo^^^^^^V^il^p^oac^dr^s^ot^hMac- 
terizlT robustness mdZts nSy ^extend to msmethoS 

The u method can be immediately used to compute aeroservoelastic stability ” ar ®' ns „ 

ju^ll co^utes^a wors^ase stability^a^m^^^nmewo^allows flutteAnd aeroservoelastic stabil- 
ity margins to be easily computed from the same model. 

9.4 Nonlinear Limit-Cycle Oscillation Analysis 

The flutter phenomenon analyzed in *2^ 
operators. Another phenomenon known as a y , structura i dynamics, such as nonlin- 

*e ae ^ namiCS ^ ^ SeParatl ° n ^ 

shock movement. ’ 

syste^s^L^ear°o[wr^re^a^b^used^o^l^un^mK^^n^^^ ^^1 amounts of^mrnlinMnty. 
rSr^^^^bec^ these phenomena, which can he 
destructive to the aircraft, are actually stable in the sense of Lyapunov. 

An »■ the^linear m"at“ 

~ » used to generate a nonhnear operator to assoc, ate with the .inear 
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operator through feedback that allows the closed-loop system to model the observed nonlinear behavior 
of the aircraft. Uncertainty descriptions are also generated for the linear and nonlinear operators, and 
robust stability theory is used to compute a worst-case stability margin. 

'7'} ■T'l 

Two approaches are investigated to analyze a nonlinear aeroelastic system. ’ The first approach 
uses multiresolution wavelets to identify nonlinear behavior in the flight data . 74 These wavelets present 
time and frequency information that allow a more accurate representation of the observed dynamics . 75 
Another approach is based on nonlinear system identification in the LFT framework . 76 This approach 
will simultaneously generate a nonlinear operator and an uncertainty description for use with robust 
stability analysis. 
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CHAPTER 10 
CONCLUSION 

A structured singular value (p) analysis method of computing flutter margins has been introduced. 
This method analyzes robust stability of a linear aeroelastic model with uncertainty operators. Flight data 
can be used to formulate the uncertainty operators to accurately account for errors m the model and the 
observed range of aircraft dynamics caused by time-varying aircraft parameters, nonlineanties, and flight 
anomalies such as test nonrepeatibility. The p-based approach computes flutter margins that are robust or 
worst-case, with respect to the modeling uncertainty. 

Nominal and robust flutter margins have been computed for the F/A-18 Systems Research Aircraft 
(SRA) using p and p-k methods. The similarity of the nominal flutter margins demonstrates the p method 
is a valid tool for computing flutter instability points and is computationally advantageous. Robust flutter 
margins have been generated with respect to an uncertainty set generated by analysis of extensive flight 
data. These margins are accepted with a great deal more confidence than previous estimates because 
these margins directly account for modeling uncertainty in the analysis process. The robust flutter mar- 
gins indicate the desired F/A-18 SRA flight envelope should be safe from aeroelastic flutter instabilities; 
however, the flutter margins may lie noticeably closer to the flight envelope than previously estimated. 

This method introduces p as a flutter margin parameter that presents several advantages over tracking 
damping trends as a measure of a tendency to instability from available flight data. The p extends to 
provide information about model sensitivity and is equally suitable for aeroelastic and aeroservoelastic 
analysis. The predictive nature of p makes it particularly attractive for flight test programs because a 
robust stability margin can be computed using data from a stable test point. 

Dryden Flight Research Center 

National Aeronautics and Space Administration 

Edwards, California, October 20, 1997 
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APPENDIX A 

UPPER BOUND FOR THE STRUCTURED SINGULAR VALUE 


The structured singular value, p, is a difficult quantity to compute and appears to be a nonpolynomial- 
hard (NP-hard) problem. 7 Upper and lower bounds have been derived that can easily be computed. The 
lower bound can be solved by a power iteration that extracts the smallest destabilizing perturbation. 20 The 
upper bound can be posed as an optimization problem to present a conservative bound on the worst-case 
stability properties of the system. This paper is primarily concerned with computing a worst-case stability 
analysis, so this appendix concentrates on the optimization formulation to compute a p upper bound. 

Question A.0.1: What is an upper bound for p with respect to a general structured uncertainty 
description that can be computed? 

A simple upper bound for |i is immediately obvious from definition 3.3.1. 


p(P)<c(P) (Al) 

This upper bound is essentially the robustness as measured by the small gain theorem and can be highly 
conservative when the uncertainty operator is structured. The conservatism in this upper bound can be 
reduced by considering two properties of singular values. 

Property A.0.2: The maximum singular value has a multiplicative property. 

a(P)c(P _1 ) = 1 (A2) 

Property A.0.3: The maximum singular value has a triangle-inequality property. 

ct(DP) < a(D)a(P) (A3) 

These properties demonstrate the maximum singular value of a scaled matrix is less than or equal to 
that of an unsealed matrix. 

Lemma A.0.4: Given matrix Pe C nXn and invertible D e C n * ", then 

a(DPD ')<a(P) (A4) 

Recall the definition of the uncertainty structure given in section 3.2. 


(A diag ( 5 1 l R } ’ 5 m r R m ’ 5 1 ^ A l’ •••> A p ) 


,R 


: Sj e R, 5j e C, Aj e C ' 


C: x C; 


(A5) 


A set of scaling matrices is defined with structure similar to the uncertainty structure. The set Z> 
contains blocks to scale with every complex and real uncertainty block. 20 


81 


V = |d = diag Dm' D ?’ •"’ D n ,<! 1 'c,' •••’ d P I «p) 

: D f s C R ' X C c ' X Cl , d, d C, D = D* > o| 
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The set of scalings $ affect only the real parametric uncertainty blocks. 

f R: xR j _J 

£ = G = diag (Gj, .... G m , 0, .... 0) : Gj 6 C , G = G 


(A6) 


(A7) 


Several advantageous properties of these scalings exist. 

Property A.0.5: Symmetric factors and inverses of scalings are valid scalings. 


_! 1 

D -1 , D 2 , D 2 e T> for every D e V 


(A8) 


Property A.0.6: Every scaling DeVandGe# commutes with every uncertainty A e A. 


DA = AD 
GA = AG 


(A9) 


The Li is shown to be invariant under transformations using the scaling matrices V. 20 This invariance 
arises as a direct consequence of Properties A.0.5 and A.0.6, which demonstrate A = VAD . Thus com- 
puting \i with respect to A is equivalent to computing \l with respect to DAD or a G • 


upper bound for (I is formulated using these scaling matrices. 


li(P)< inf o(DPD ^^(P) 

DeZ> 


(A 10) 


This upper bound is less conservative than the unsealed upper bound because of the final mequah ty 
that is a consequence of lemma A.0.4. This upper bound uses a search over all valid scaling m ices ^to 
find the optimal scaling matrix that reduces the conservatism in the robustness computation. This upper 
bound can be formulated as a maximum eigenvalue condition as demonstrated in lemma A.0.7. 


Lemma A.0.7: The constant matrix P is robustly stable and p(P) < 1 with respect to the set A if 


or equivalently 


inf ||D 2 PD 2 IU < 1 

DeZ> 


inf )i(P*DP-D)<0 

DeZ> 


(All) 


(A 12) 


This upper bound allows scalar uncertainty parameters to lie within a norm-bounded disk centered in 
the complex plane that can be conservative if these parameters are real. The conservatism can be reduced 
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by introducing the second set of scaling matrices, * that affect only the real parametric uncertainty 

blocks. These scalings restrict the uncertainty set to lie within a norm-bounded disk covering the real axis 

b " and T 1,but that allowed t0 be centered off the origin. A new upper bound is formulated 
with both Z> and ^ scalings. /8, /y 


p(P) < Jmax (0, a*) 

a* = inf min {a : P*DP + y'(GP - P*G) - aD < 0} 

DeV aeR 

G€^ (A13) 


This upper bound is presented as an optimization in lemma A.0.8. 

Lemma A.0.8: The constant matrix P is robustly stable, and |l(P) < 1 with respect to the set A, if 


inf {a : ^(P*DP + j(GP- P*G) - aD) <0} < 1 

DeT> 

Ge<? 


(A 14) 


The ^matrices are zero if no real parameter uncertainty blocks exist in the system, so the condition in 
emma A.0.8 reverts to the condition in lemma A.0.7. When the ^ matrices are nonzero, the additional 
terms iyhe new upper-bound condition arise by using an additional constraint on the signals in the 
system. Consider a signal in the closed-loop system relating an uncertainty signal. The system 

equation takes on the following form if the system is well-posed, ignoring the extraneous error and 
disturbance signals. 


APx = x 


(A 15) 

Allow the uncertainty to be a single repeated real parameter, A = 8 1 with 8 e R, and take the complex 
conjugate transpose of the entire equation. ^ 


8x*P* = x* 

Combining this formula with a scaling G e # leads to the following. 

x *GPx = 8x*P*GPx = x*P*G(8Px) = x*P*Gx 
The following matrix is derived from the right and left sides of this equation. 

x*(GP- P*G)x = 0 

Consider the matrix term of the complex uncertainty robustness upper bound scaled by x. 

x*(P*DP - aD)x < 0 

Combining the two equations leads to the new upper-bound condition. 

x*(P*DP + j(GP - P*G) - aD)x < 0 


(A 16) 


(A 17) 


(A18) 


(A 19) 


(A20) 


83 


This p upper hound is directly re.aied ,0 robust 

rent™ The t> scaling matrices use structure in ihe u ^ ^ ustness condit ion is also related to a 

The p upper-bound conditions * n ^^^^^^ghouT^h^derivahon^In practice^the 

^ tr; v ^r t"e 

case upper bound. 

. . . • P / ; rn \ p is robustly stable and |i(P) < 1 with 

Lemma A.0.9: The transfer-function matrix P (j co) e <*oo 

respect to the set A if 


max inf {a : l(?U «)*DPO •) + MW “> - P * ( ^ )G) ' « D > S 01 < 1 

(D D e t> 

Ge # 


(A21) 


* nnnpr hound bv finding optimal scaling matrices 

Established commute with linear time-invariant (LTI) 

at distinct frequency points. Frequency vary S * modeling errors. Restricting the scaling 
uncertainty that may not accurately cap ure „ consideration of uncertainty operators with 

matrices ,o he cor is a necessary and stiff, • 
linear time-varying ( , with com plex LTV uncertainty for any block struc 

cient condition for robust stability of p m§ ^ ith comp i e x LTV uncertainty has been 

ture. 8 Computing robustness measures fo y , 83 R lated papers have studies of robust- 

addressed for finite-energy signals 82 and for fimte-gam signals^ 84 

ness conditions for systems with some type of rate bounds on the LTV uncertainty. 

•ii • ~ Knth T TT and LTV uncertainties blocks to accurately describe the 
modelin^errors. ZZ7'c —ties are associt^ wi* 

tured uncertainty operators. 


A xv (t){ Aft) = diag(5?I R| 6 P I Rm . 5^I Ci 8^c n ' A A P 

~ C i XC i] 

•.sfcOeReSfcOeCVOeC 


(A22) 


(C0)|a(C 0) = diag^I Ri , •••> 5 m I R m ’ 5 l I C,’ •••’ 8 n I C„’ A l’ •” ,A P 

r c i xc i] 

:5f(to)eR,8i(a)eC,A i (CD)€C 


(A23) 


m „ , TV nnrprtaintv oDerators that are allowed to change infinitely fast as a 

,f m T e V 0 ,™chrhave invesiafed limits on the rate variation of complex time-varying 
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parameters that are not addressed in this paper. 84 Correspondingly, the set A TI contains all LTI uncer- 
tainty operators that are allowed to vary with frequency. These sets are allowed to enter the problem as 
shown in figure A.l. 



Figure A. 1 . Block diagram with LTI and LTV uncertainty. 


Associated with these uncertainty descriptions are scaling matrices. Constant scaling matrices, Z> TV , 
are structured like T> to associate with the time -varying uncertainty, A TV . Frequency-dependent scalings, 
Z?Ti (©). are associated with the time-invariant uncertainty, A TI . The block structure of each of these 
sets matches the block structure of the corresponding uncertainty set. 

The uncertainty and scaling matrix sets can be combined into single sets for ease of notation. 

A = diag{A TI (to), A xv (t)} 

“D — diag{D TI (co), D-py } (A24) 

$ — diag{G TI (co), G-py } 

Lemma A.0.10 describes the robustness upper-bound condition. 

Lemma A.0.10: Given some a > 0 and norm-bounded uncertainty A such that 

max a(A TT (co)) < - max a(A TV (t)) < - 

“ TI a t TV a (A25) 

then the system described by transfer function P (j to) is robustly stable with respect to A if scaling matrices 
D e Z> and G e # exist such that : 

max X(P(yG))*D(a>)P(y'co) + j(G(y'co)P(y'co) - P(yco)*G(oo)) - a 2 D(co)) < 0 (A26) 

The search for optimal scaling matrices that achieve the lowest robustness value possible can be for- 
mulated as an optimization to minimize a subject to a constraint function associated with an eigenvalue 
condition. An additional constraint arises from restricting the Z> scaling matrices to be positive definite. 
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min a 

CO 


subject to -M*DM - j (GM - M*G) + a D>0 

J D > 0 


(A27) 


This upper-bound condition is written in this format to ease understanding of the function an 
develop solution algorithms. Specifically, the minimization problem as shown above with the eigenvalue 
constraints is a linear matrix inequality (LMI). The LMI is a convenient framework to use because com- 
mon properties global to all LMI problems exist. 85 The mam property states an LMI is a convexfunct 
of the viLbles that allows solutions to be computed using standard convex optimization algorithms The 
ellipsoid method and the interior point algorithm of the method of centers are shown to be efficient for 

LMI solutions. 86-88 

This optimization involves a search over frequency. Only a finite set of frequency points are used in 
practice, and engineering judgment must be used in choosing the set of frequency points. Testing he 
robustness condition at more frequency points will increase the accuracy of theupper bound, bu e 
number of free parameters will rise accordingly and increase computational time. 
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APPENDIX B 

NOMINAL AEROELASTIC MODELS IN THE STRUCTURED SINGULAR VALUE 
FRAMEWORK WITH KNOWN AERODYNAMIC STIFFNESS AND DAMPING 


This paper presents an aeroelastic model formulated by approximating the unsteady aerodynamics as 
a state-space system. The unsteady aerodynamic forces are treated as an unknown system with no 
specialized structure. This appendix considers formulating a state-space system to the unsteady aerody- 
namics with a specified structure. The resulting aeroelastic model is parameterized around dynamic pres- 
sure, and the structured elements of the aerodynamic forces appear in the model. 

The general formulation shown in section 4.1 approximates the unsteady aerodynamic forces as a 
state-space system. This unstructured formulation does not explicitly include additional information 
regarding the aerodynamics that can be obtained from sources other than the computational model. In par- 
ticular, the steady aerodynamics and first-order derivatives are often measured from wind-tunnel testing. 
Two matrix terms are introduced to the unsteady aerodynamic formulation to include this information. 


Q(s) = Aq + sA, + 


A Q Bq 
C Q Dq_ 


- 1 , 


= A 0 + sAj + Dq + Cq(sI - Aq) B 


(Bl) 


2 This system is quite similar to Karpel’s form. 29 The main difference is the absence of an additional 
s A 2 term. This term is typically not available from experimental data and is not explicitly represented in 
this formulation. 


The state-space elements, {Aq, Bq, Cq, Dq}, are computed using standard system identification 
algorithms. The known matrix elements are subtracted from the unsteady aerodynamic force data, and the 
resulting system, Q, s (s), is approximated by the state-space elements. 


Q; S (S) = Q(s)- Aq-SAj 




(B2) 


Given the number of generalized states, n, such that A Q , Aj e R 


11 /v 11 


states, nQ, define Aq e R 


n O xn o n n xn nxri,. 

, Bq6R Q , C Q e R Q ,andD ( 


and the number of aerodynamic 
„n x n 

e R as the elements of the 


state-space system approximating Q ss (s). 


The aeroelastic model is extended to include the additional Aq and Aj terms. Consider the signal, y, 
generated by an input to the state-space portion of force matrix. Define the signal, x e r" q , as the vector 
of aerodynamic states. 


y = q s (s)ri 


<=> 

X 


A Q Bq 


X 


y 


C Q D Q_ 




(B3) 
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Formulate the aeroelastic differential equation using x. 


0 = Mf\ + Cr\ + Kti + qQ(s)r\ 

= Mfi + Cf\ + Kr| + q(A 0 + sAj + Dq + C Q (si - A q ) B Q )T[ 
- Mf\ + (C + q Aj )fl + (K + <?A 0 )ri + ?Qss (s) 1 ! 

= Mri + (C + + ( K + 4A 0 )t| + W 

= Mti + (C + q A^fi + (K + ^A 0 )r| + q(C q x + D q ti) 

= Mil + (C + ^Aj)!! + (K + qA Q + ^D q )r\ + qC Q x 


n a^c sta^^ following *"* ' differential eqUati ° nS - 


formulated using the generalized states, ti and Tl , and the unsteady aerody- 

. . .1 r 11 : « /4«ffar , ontlOl #*ni 1 OT"l 


11 


o I 0 

-M'Vk + ^Aq + ^Dq) -M'hc + qA,) -qM C Q 

11 

ii 


Tl 

X 


b q 0 A Q . 

X 


This stale-space aeroelastic can be parameterized around dynamic pressure. Consider an additive 
perturbation to dynamic pressure, 8- e R. 


4 = Vo +b q 


(B6) 


Separate terms in the system dynamics that involve 8- . 


0 = Mfj + (C + ^Aj)"^ + (K + <?A 0 + «D q )T| + qCqX 

= Mii + [(C + ^A ] )f\ + (K + ^A 0 + gD Q )Ti + 9 o C Q x] 

+ S- [A.fl + AqT| + Dq T l + Cq x] 

= f^ + [M -1 (C + § o A 1 )fi + M" 1 (K + ^ o A 0 + ^,D Q )Ti + ^M C Q x] (B7) 

+ 8- [M _1 A,fl + M _1 A 0 r| + M 'DqTI+M C q x] 

= f\ + [M -1 (C + A|)il + M -l (K + A q + ^ 0 D q )ti + C q x] + S^z 

= fi + [M _1 (C + ^A 1 )f|+M“ 1 (K + ^A 0 + ^D Q )tl + ^M C Q x] + w 


The signals z and w are introduced into this formulation to associate the perturbation in dynamic 
pressure to the nominal dynamics in a feedback manner. The signal z can be generated as an output 
plant because it is a linear combination of states. 

(B8) 


z = M 1 Ajfj + M AqT) + M D ( 


,r| + M ’C ( 
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The signal w is related to z by the dynamic pressure perturbation. 


The state-space aeroelastic model for nominal stability analysis in the p framework is formulated 
using the state-update matrix. The matrix is determined by the dynamics at the nominal dynamic pressure 
and the additional input and output signals used to introduce perturbations to the dynamic pressure. That 
perturbatron, 5-, is not an explicit parameter in the state-space model because the perturbation only 
affects the plant through a feedback relationship as determined by the signals z and w. Define the transfer 
function P(s) generated by state-space matrices such that z = P(s)w. 




0 

I 

0 

0 

Ti 


-M ’(K + ^Aq + ^Dq) 

-M _1 (C + ^A,) 

-m l c Q 

-I 

X 


b q 

0 


0 

z 


^ ( Aq + Dq ) 

M _1 Aj 

M~' C q 

0 


(BIO) 


Figure 4.1 shows the feedback interconnection between the perturbation in dynamic pressure and the 
nominal plant model parameterized around that perturbation. This interconnection is a linear fractional 
transformation, and the small gain condition of lemma 2.4. 1 or the structured singular value condition of 

_ C ^ em 3 ' 3 ' 3 can j be directl y applied to analyze stability with respect to a variation in the flight condition 
q. 1 his system reduces to the parameterized model presented in section 4.2 if A Q = Aj =0. 
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